Search for Ferromagnetism in doped semiconductors in the absence of transition 

metal ions 
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In contrast to semiconductors doped with transition metal magnetic elements {e.g. Gai-^Mn^As), 
which become ferromagnetic at temperatures below ~ 10 2 K, semiconductors doped with non- 
magnetic ions {e.g. silicon doped with phosphorous) have not shown evidence of ferromagnetism 
down to millikelvin temperatures. This is despite the fact that for low densities the system is 
expected to be well modeled by the Hubbard model, which is predicted to have a ferromagnetic 
ground state at T = on 2- or 3-dimensional bipartite lattices in the limit of strong correlation 
near half-filling. We examine the impurity band formed by hydrogenic centers in semiconductors at 
low densities, and show that it is described by a generalized Hubbard model which has, in addition 
to strong electron-electron interaction and disorder, an intrinsic electron-hole asymmetry. With 
the help of mean field methods as well as exact diagonalization of clusters around half filling, we 
can establish the existence of a ferromagnetic ground state, at least on the nanoscale, which is 
more robust than that found in the standard Hubbard model. This ferromagnetism is most clearly 
seen in a regime inaccessible to bulk systems, but attainable in quantum dots and 2D heterostruc- 
tures. We present extensive numerical results for small systems that demonstrate the occurrence of 
high-spin ground states in both periodic and positionally disordered 2D systems. We consider how 
properties of real doped semiconductors, such as positional disorder and electron-hole asymmetry, 
affect the ground state spin of small 2D systems. We also discuss the relationship between this 
work and diluted magnetic semiconductors, such as Gai_ I Mn 1 As, which though disordered, show 
ferromagnetism at relatively high temperatures. 



I. INTRODUCTION 



Originally proposed in the early l^O o 1 ' 2 ' 3 ^, the 
Hubbard model combines tight binding hopping be- 
tween nearest neighbors on a lattice with an on-site 
Coulomb repulsion between electrons in the same or- 
bital state. Though it is one of the simplest inter- 
acting models, its on-site intra-orbital correlations are 
believed to be the most important source of correla- 
tions in solids. Indeed, the Hubbard model displays 
great diversity of transport and magnetic properties, 
giving rise to insulating, metallic, and superconducting 
phases as well as ferromagnetic (FM), antiferromagnetic 
(AF) and paramagnetic spin order. It has been used 
to study a wide range of correlated systems, including 
Mott-insulator oxides' high-T c superconductors'^! 9 - 
organic materials ) 10 ' 11 ! 12 -\/3-adlayer structures,— vana- 
dium oxides ; 14 ' 15 nickel sulphide-selenide alloys ) 16 ' 17 ' 18 
hydrogenic centers in doped semiconductor s 19 ' 20 , and 
quantum dots^i Such great interest and applications 
have resulted in analyses of the model on different 
lattices ) 22 ' 23 with multiple^ 4 - and degenerat o 25 ' 26 bands, 
and with binary alloy disorder. — Many studies restrict 
themselves to the infinite U/t limit ) 28 ' 29 which can be 
realized most effectively in optical lattices' - but can be 
approached in semiconductor systems as well. We will be 
concerned with the case of semiconductors doped with 
shallow hydrogenic impurities. Here the model is partic- 
ularly appropriate at low densities {i.e. in the insulating 
phase, where carriers are bound to a few sites and the 



Coulomb interaction is large compared to the kinetic en- 
ergy). In this low density limit each site is treated as an 
effective hydrogen atom with a corresponding effective 
Rydberg and Bohr radius: 
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where m* is the effective mass in the appropriate band 
and e is the dielectric constant of the host material. 
In doped semiconductors, typically e ~ 10 — 20 and 
to* is 0.05 to 0.5 times the free electron mass, so that 
a*, ~ 10 - 500 A and Ry* ps 1 - 50meV. Since the Ry* is 
usually much smaller than the bandgap of the host semi- 
conductor, the lattice lacks low-energy electronic excita- 
tions on the energy scale of the impurity electrons and 
essentially plays the role of an inert vacuum. Realistic 
effects like valley degeneracy and mass anisotropy must 
be included for quantitative calculations but are unneces- 
sary for the qualitative phenomena of interest to us. 31 ' 32 
We will assume that all relevant energy scales are much 
smaller than the gap between the lowest and higher or- 
bital states on an isolated dopant, so that we need only 
care about the Is orbital of each dopant, which consists of 
two electronic spin-degenerate states at energy denoted 
Eq. A hydrogenic center, like a hydrogen atom, is known 
to bind up to two electrons! 3 - 3 - With a single electron the 
problem is that of atomic hydrogen {H), and the electron 
is bound with 1 Ry* . The two electron case corresponds 
to the H~ ion, which has a spin singlet ground state 
bound by 0.0555 Ry* 

We begin with a review of the Hubbard model and 
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its properties on a lattice in section [TTJ The absence of 
certain magnetic properties, namely ferromagnetism, in 
real materials leads to a discussion of disorder and re- 
veals the need to incorporate it into the model. This is 
done in section ITTT1 where we motivate and define a model 
appropriate for doped semiconductors. The parameter 
ranges of interest for this model are also given in section 
IIII1 along with details of the model's solution. Results 
on finite lattices, selected symmetric clusters, and small 
random clusters are presented in section IIVI Large sys- 
tems of random impurities are treated in section [V] by 
dividing them into smaller clusters which can be solved 
exactly. Section IVT1 highlights our major conclusions and 
discusses topics for continued work. 



II. BACKGROUND: THE HUBBARD MODEL 



A. Definition and general properties 
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FIG. 1: (Color online) Schematic figure showing a system at 
half-filling (a), and slightly above half-filling (b). At half- 
filling the lower impurity band is completely full and there 
is a gap to charge excitations. Above half-filling there are 
electrons present in the upper (unfilled) band that can act 
as carriers if they occupy extended states (as they do in a 
lattice). Note also that each band's density of states N(E) is 
not actually semicircular, but drawn this way for convenience. 



The Hamiltonian of the Hubbard model on a lattice 
with N s sites is given by: 

n =- t J2 ( C *W + h.c) + U ]T n^n iL (2) 

{i,j)a i=l 

where i and j range from 1 to N s , and the first sum is 
over all distinct nearest neighbor pairs. Operators c\ a 
and eta- create and annihilate, respectively, an electron of 
spin <7 G {T: 1} on site i, and satisfy canonical fermion 
anticommutation relations 

{ct^c^ j = fe<W (3) 

{4.4.'} = o (4) 

{c^c^/j = (5) 

for any i = 1 . . . N s and a G {Til}- The number operator 
rii a = cl a Ci a , and has eigenvalues and 1. The Hilbcrt 
space of each site has as a basis the four states where 
the site is occupied by an up-spin, a down-spin, neither, 
or both. The total Hilbert space is the direct product 
of site Hilbert spaces, and therefore has dimension 4^° . 
The parameter t is the quantum mechanical hopping am- 
plitude between (nearest-neighbor) sites, and U is the 
strength of the on-site Coulomb repulsion. We include a 
minus sign in front of the kinetic term, so that for the 
familiar example of the tight-binding model with hydro- 
genic wavefunctionspS, t(r) = 2(1 + r/as) exp(— r/ae) is 
positive, and restrict ourselves to U > (repulsive inter- 
action) . Note that the eigenstates of each term indepen- 
dently are trivial: they are states of definite momentum 
when U = and states of definite position when t — 0. 
Thus, inherent in the Hubbard model is a competition 
between the extended (wave-like) and localized (particle- 
like) nature of the electrons, and there is no clear classical 
analogue. 



We will be primarily concerned with the strong corre- 
lation limit U S> t, so that at half-filling (i.e. one electron 
per site) the single particle (charge) spectrum has a gap, 
and the system is insulating (Fig. [Ua)). Nonetheless, 
the system can have low lying spin excitations. At large 
U/t and half-filling, the Hubbard model at low energies is 
effectively a Heisenberg model^l in which fermionic elec- 
tron operators are represented by spin operators. On a 
bipartite lattice, the Heisenberg model, given by Hamil- 
tonian Ti-Heis = JJ2<ij> Si ' Sji nas an AF ground 
state. 38 Away from half-filling, where there are carriers 
(see Fig. [IJb)), one must use a more general low-energy 
theory that includes a kinetic term, the t — J mode h 39 i 40 

U t j = -t y~] nl - n i& )c\ a c ja (l - rijff) +h.c.^ 

<ij>o 

+ J ~ \ n i n i) ■ 

<ij> ^ ' 

Note that the Hamiltonian operates on the restricted 
Hilbert space which excludes doubly-occupied sites. 
From the Heisenberg and t — J models we see that 
the inclusion of electron-electron interactions results in 
an AF exchange interaction ~ Si ■ Sj, where Si — 

c \ a a ap c iP- The exchange term is due to the fact 
that virtual hopping of electrons between neighboring 
sites is allowed when their spins are oppositely oriented 
but not when their spins are parallel (as in a FM config- 
uration), as shown in Fig.[2^i 

The Hubbard Hamiltonian can also be written in 
terms of spin operators using the identity J2i (>Si) = 
Si (l n i - I'M 71 *!)' where n t = n^+nn, casting Eq. © 
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FIG. 2: (Color online) Diagrams which graphically tell the 
origin of the AF exchange interaction term of the t — J model: 
in an AF configuration (left) electrons and virtually hop to a 
neighboring site and back (shown by the arrows), resulting in 
a net lowering of the energy by second order perturbation the- 
ory. In a FM configuration (right), however, Pauli exclusion 
forbids such virtual processes, and the system cannot lower 
it's energy in this way. 



into the form: 

H = ~ f E (4r<* 



h.c. - 
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(7) 



where N e is the total number of electrons. This form 
clearly shows the total spin SU(2) invariance of the Hub- 
bard model, and also that when U > the interaction 
energy is lowest when the total spin on each site is max- 
imized, suggesting the existence of ground states with 
high values of total spin at large U. On a bipartite lat- 
tice with disjoint sublattices A and B, the sign of t can 
be changed via the transform: 



-Ci a if i G A 
-Ci a if i € B 



(8) 



which does not change the canonical commutation rela- 
tions and thus leaves the spectrum invariant. The Hub- 
bard model also possesses particle-hole symmetry on a 
bipartite lattice, where U maps to — U and total charge 
is interchanged with total spin (for a detailed explanation 
of symmetries in the Hubbard model, see Ref . l42l) . 

Even when it is applied to simple systems {e.g. 1-, 2-, 
and 3-dimensional lattices) , the Hubbard model yields in- 
teresting and non-trivial properties, seen through the na- 
ture of its excitations, density of states, spectral weight, 
transport, and optical and magnetic behavior i 43 i 44 ' 45 i 46 
Here we will concentrate on the nature of magnetic cor- 
relations in the ground state, which are then used to con- 
struct the ground state (i.e. T = 0) phase diagram. 



B. Magnetic Properties 

The magnetic properties of Hubbard systems can be 
very rich due to competition between two or more mag- 
netic phases. Consider the Hubbard model at half-filling 
on a bipartite lattice, where there is no classical magnetic 



frustration, and let U/t be large. The model's quantum 
ground state is a superposition of "Neel antiferromagnet 
states" where spins on each sublattice are aligned and op- 
positely oriented to those of the other sublattice as well 
as "spin-flip states" which differ from the Neel AF states 
by exchanging one or more pairs of spins (states with a 
greater number of flips occur with lower weight). In other 
words, the ground state is a superposition of states with 
long-range Neel order. Since U is large, the t—J approxi- 
mation (Eq. ©) is valid, introducing an exchange energy 
J ~ t 2 /U between neighboring spins. The kinetic term 
of the t — J Hamiltonian does not play a role since at half- 
filling there are no mobile carriers. Thus, at half-filling 
the exchange directly gives rise to an antiferromagnet. 
When the system is above or below half-filling, however, 
the kinetic term plays a competing role by favoring a 
ferromagnetic spin configuration. This is so because as 
carriers hop from site to site they do not disturb an un- 
derlying FM spin configuration, whereas they necessarily 
scramble an AF one (see Fig. [3J . This scrambling leads 
to an unfavorable increase in energy, and thus the pref- 
erence for ferromagnetismi 47 ! 48 Relative to an AF state, 
a FM system with carrier (electron or hole) density S 
gains kinetic energy of order td due to carrier dereal- 
ization and loses magnetic energy of order is J = 4t 2 /U. 
Thus, at a fixed small 8, when U is large enough, tS ^> J, 
and the system prefers a FM configuration over the AF 
one because it allows carriers to be less confined. Under- 
standing the applicability and validity of this argument, 
and more generally the factors that govern the magnetic 
competition found in the Hubbard model, has been the 
topic of much work. Indeed, it has led to most (if not 
all) of the rigorous results that are known concerning the 
Hubbard model. 

Despite the apparent simplicity of the Hubbard model 
and voluminous literature surrounding it, few rigorous 
theoretical results have been proven about it. Most 
striking among them is the result of Nagaoka, 49 which 
states that in the infinite correlation limit U/t — > 00, 
the Hubbard model on certain finite lattices of dimen- 
sion d > 2 with periodic boundary conditions, t < 0, and 
a single hole (away from half-filling) , has a FM ground 

state (i.e. the total spin S 2 , where S = J^i^i, attains 
it's maximal value). This result, dubbed the Nagaoka 
Theorem, applies to most standard lattices, including 
the square, simple cubic, triangular, kagome, bec, and 
fee (hep) i 49 i 50 In the case of bipartite lattices, such as 
the square, simple cubic, and bec, t can be taken pos- 
itive (the physical sign in the tight binding model) by 
the transform of Eq. ([5]). This can be understood from 
the preceding discussion of a bipartite system, where, 
upon setting U — 00, the criterion t5 3> J is satisfied 
for any S > and thus only a single hole is needed 
to produce a FM ground state. Even though this cri- 
terion also predicts ferromagnetism for a finite density 
of carriers at large U, a rigorous result even for the case 
of a few holes has proved difficultj 51 ' 52 Along with the 
rigorous proofs in Nagaoka's and Thouless' work, sim- 
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FIG. 3: (Color online) Diagrams showing why the kinetic term 
favors a ferromagnetic state: in a) the down spin on the single 
doubly-occupied site can move freely without disturbing the 
underlying FM background. However, if the background is 
AF as in b), motion of electrons on doubly-occupied sites 
scramble the Neel order. Diagram c) shows the result of the 
doubly occupied site in b) moving two sites to the right. 



pier and more modern mathematical proofs are given by 
Tian 53 and Tasakii^. Another rigorous theorem regard- 
ing magnetism in the Hubbard model by Lieb^ states 
that a half-filled bipartite system whose sublattices have 
different numbers of sites will have an unsaturated FM 
ground state. It later became clear that the tendency to- 
ward ferromagnetism was due to the single particle den- 
sity of states being dispersionless, or flat, at the center of 
the band (the Fermi level at half- filling) . Later, results 
of Mielke^ and TasakiSi generalized this idea to charac- 
terize a broader class of half-filled systems with disper- 
sionless single-particle spectra and saturated FM ground 
states, said to exhibit "flat-band ferromagnetism." We 
hasten to point out, however, that half-filled Hubbard 
systems are generally antiferromagnetic (when on a bi- 
partite lattice) or paramagnetic, and that completely or 
nearly flat bands should be viewed at least as a non- 
generic case. This fact underscores the surprising result 
of Nagaoka, which describes the transition from an an- 
tiferromagnet to a ferromagnet upon the addition of a 
single hole or electron. 

The bulk of this section investigates the possibility of 
saturated "Nagaoka ferromagnetism", and it is worth- 
while at this point to consider the progress of past 
work toward understanding the phenomenon. The topic 
has generated sizable interest, since the Nagaoka The- 
orem is at the same time striking and of only lim- 
ited use, saying nothing about the thermodynamic limit 
where there is relevance to experiment. Many theoretical 
studie o 28 ' 29 i 58 ' 59 i 60 work in the large or infinite U limit, 



which is where saturated FM is most likely to occur. In 
the U = oo limit doubly-occupied states are eliminated 
from the Hilbert space, which then has a dimension that 
scales as 3^ - substantially less than 4^ and thereby 
a great relief for numericists! Indeed, much computa- 
tional work has been done setting U = oo, including one 
which relates the system far below half-filling to one of 
hard-core bosons^ 9 - 

Investigating the existence, extent, and stability of the 
Nagaoka state has established several conditions that are 
known to favor a stable, saturated FM ground state. By 
considering the stability of the fully polarized state to 
a single spin flip, it is show n 22 i 61 i 62 that an asymmetric 
density of states with a peak at the appropriate band 
edge (lower edge of the upper Hubbard band if doping 
above half-filling; upper edge of lower Hubbard band if 
doping below) is one such condition. This makes intu- 
itive sense, since having large density of states at the 
Fermi level diminishes the kinetic cost of filling additional 
single-particle electronic states and causes the large U, 
which favors spin alignment, to prevail. This generalizes 
the condition of a flat band discussed earlier, in which 
the density of states is infinite at the band edge. From 
the geometrical optimization of finite Hubbard clusters, 
Pastor et alj 61 i 63 find that saturated ferromagnetism co- 
incides with clusters which are non-bipartite and have a 
large number of frustrated "tight" triangular loops. They 
also find that doping clusters above rather than below 
half-filling yields a density of states with higher weight 
at the band edge, and leads to FM ground states. The 
asymmetry with respect to doping and correspondence 
between magnetism and triangular loops is corroborated 
by our results, and appears to be a quite general feature 
with important experimental ramifications, which we dis- 
cuss in more detail in section TlV CI It is also known that 
adding back into the single-band Hubbard model physical 
interactions that it neglects, particularly a (direct) fer- 
romagnetic Heisenberg exchange interaction, can be im- 
portant for stabilizing ferromagnetism near half-filling for 
finite t/i 64 ' 65 i 66 Additionally, the next-nearest-neighbor 
(NNN) hopping amplitude t' is believed to play an im- 
portant role: decreasing t'/t (especially below zero) sta- 
bilizes saturated FM to higher hole-doping in the U = oo 
Hubbard model on a square latticed In one-dimension, 
where the Lieb-Mattis theorem^ 8 - forbids FM in the stan- 
dard Hubbard model, the addition of (NNN) hopping t' 
such that t'/t < (t > 0) results in a widespread FM 
phased Alternatively, the generalization to a multi-band 
model (appropriate for many transition metals) with a 
ferromagnetic exchange interaction between electrons in 
different orbitals ("Hund's rule couplings") also abets 
the stability of a FM stated (The inclusion of multiple 
bands, however, is not crucial to FM stability in 2 and 3 
dimensions;—) We do not consider either of these routes, 
and restrict our study to a nearest-neighbor model with 
one orbital and on-site Coulomb interaction. 

It is important to remember, however, that saturated 
ferromagnetism in the Hubbard-like models is not ubiq- 
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uitous. Other work has shown it to be a subtle effect, 
depending on dimension and lattice geometry. For in- 
stance, another rigorous result, due to Lieb and Mattis,— 
proves that in finite one-dimensional systems with zero- 
wavefunction or zero-derivative boundary conditions, the 
ground state must be a singlet (no spin polarization). 
More recently Haerter and Shastry 72 have shown that on 
the frustrated triangular lattice an itinerant hole actu- 
ally helps to produce an antiferromagnetic ground state. 
They suggest that this phenomenon holds on all lattices 
with "electronic frustration," defined as those for which 
the sign of the hopping amplitude around the lattice's 
smallest closed loop is negative. (Note that Nagaoka's 
theorem only applies to un-frustrated systems.) 
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FIG. 4: Zero temperature mean-field theory phase diagram 
of the Hubbard model on a 10 x 10 square lattice (top) and 
8x8x8 (512 sites) simple cubic lattice (bottom). Doping 
(horizontal axis) is defined as the number of extra electrons 
(above half-filling) per site. 



C. Elusive Ferromagnetism 

A qualitative picture of the Hubbard model's mag- 
netic behavior at zero temperature can be obtained by 
a mean-field analysis on square and simple cubic lattices, 
which results in the phase diagrams shown in Fig.[¥] The 
phase transitions were found by comparing the ground 
state energy and spin-spin correlations of self-consistent 
mean-field calculations that were initialized in paramag- 
netic, AF, FM, and random configurations. Our analy- 
sis does not include the possibility of phase separation, 
e.g. the existence of polarons corresponding to "carrier- 
rich" ferromagnetic and "carrier-poor" antiferromagnetic 
regions. If it occurs, phase separation could substan- 
tially alter— the simple phase diagrams given here. Bar- 
bieri and Young construct phase diagrams for the large-t/ 
Hubbard model in 2 and 3 dimensions using a variational 
Gutzwiller technique^ and find phase separation occurs 
in both cases. Dagotto et al.,— however, argue based on 
their results on 10- and 16-site square lattices that phase 
separation is generally absent in the Hubbard model, at 
least at short length scales. FigureH]also agrees with the 
extensive work by Hirsch^ in two dimensions. 

We focus on the region of low-doping and large U/t 
(the top left of Fig. 2]), where there is a FM-AF transition. 
As expected, at zero doping (half- filling), the system is an 
antifcrromagnet for all values of U/t due to the effective 
exchange interaction and an absence of mobile carriers. 
As the doping is increased from zero, it is clear from 
the mean field perspective that for large enough U /t we 
expect, on some mesoscopic or macroscopic length scale, 
a transition to a FM ground state (even though its precise 
location in phase space depends on dimension as well as 
lattice structure, and requires more careful work). 

Though the stability of the Nagaoka state has been 
studied extensively and is seen to exist in the Hubbard 
model, such ferromagnetism has not been observed ex- 
perimentally. In many Mott-insulator oxides and chaco- 
genide systems this may be explained by an insufficient 
U ft to allow for ferromagnetism (and finding a natu- 
rally occurring material with large enough U/t seems un- 
likely). However, in doped semiconductors at low dopant 
densities, U ft is tunable over several orders of magni- 
tude due to the exponential dependence of the hopping 
t on the dopant spacing [e.g. t(r) ~ exp(— r/as) in the 
tight binding model] . This versatility makes doped semi- 
conductors a promising candidate for Nagaoka ferromag- 
netism, as it allows U/t to become large (~ 100 — 1000), 
achieving for all practical purposes the limit U/t — > oo 
required by Nagaoka's theorem. Despite this, the ab- 
sence of ferromagnetism in experiments on a variety of 
doped semiconductors, both uncomp ensat e d 32 ' 76 > 77 ' 78 1 79 
and compensated^ is quite clear. In these experi- 
ments, the nearest neighbor coupling, though distributed 
broadly, had a median value of 1-10K, and FM behavior 
was searched for down to much lower (mK) tempera- 
tures to probe the T = behavior. Even with the addi- 
tional hope of alternative theories that predicted ferro- 
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magnetism of Anderson localized electrons ; 81 ' 82 both un- 
compensated and compensated systems exhibited a sig- 
nificantly lower (by factors of 10-50) magnetic suscepti- 
bility compared to the high temperature paramagnetic 
Curie result, indicating that the systems were predom- 
inantly characterized by AF correlations, both at and 
slightly below half-filling. 

To understand how this can be, let us return to the re- 
quirements of Nagaoka's theorem. Placing dopants on a 
superlattice has become possible only very recently 8 ^ - in 
naturally formed doped semiconductors, including those 
of all relevant experiments, the dopants are distributed 
randomly and an important hypothesis of Nagaoka's the- 
orem is not met. Adding such positional disorder to the 
Hubbard- like description of the system turns out to be an 
important ingredient. (It is incorporated into the Hub- 
bard model by setting t — ► tij , which then depends on the 
separation nj, see section [Hj specifically Eq. ((9]) below). 

After introducing positional disorder into the Hubbard 
model, is not clear whether any of the aforementioned 
theorems and arguments for uniform systems are still 
valid (or even relevant). First, we expect a locally fluctu- 
ating carrier density, which may wash out any distinction 
between phase separation at macroscopic and mesoscopic 
length scales for such systems. Second, since the itiner- 
ancy of carriers depends on the local lattice geometry, 
when the geometry becomes spatially inhomogeneous the 
itinerancy might be suppressed, and at the very least the 
magnetic structure will show similar spatial inhomogene- 
ity. 

More precisely, it was found that the lack of low- 
tcmpcraturc fcrromagnctism in semiconductors can be 
explained by disorder localizing otherwise mobile carri- 
ers, thereby reducing the kinetic energy gain (previously 
~ t5) and destroying ferromagnetism. Bhatt and Le o 84 ' 85 
gave insight into the true nature of the half-filled (uncom- 
pensated semiconductor) case using a perturbative renor- 
malization group method tailored for the large amount of 
disorder present in the actual system. They found that 
the randomness of the dopants results in what has been 
dubbed a valence-bond glass ; 86 ' 87 random singlet^ or 
Bhatt-Lee phasei 89 ' 90 In such a state, spins pair up to 
form (spin zero) singlets (see also Ref. iQlh in a hierar- 
chical fashion, and the resulting structure and behavior 
is qualitatively different from the antiferromagnet state 
predicted on a bipartite lattice. There is no long-range 
AF order, and the magnetic susceptibility is strongly 
temperature dependent, even down to tens of millikelvin. 
That compensated semiconductors show no evidence of 
ferromagnetism 8 ^ can be attributed to the localization of 
holes on one (or a few) valence bonds, and their conse- 
quent inability to move long enough distances to disrupt 
the local magnetic arrangements. As a result, holes are 
unable to gain the kinetic energy which favors a spin- 
polarized background. Thus, even though doped semi- 
conductors give one the ability to tune U/t over several 
orders of magnitude, Nagaoka ferromagnetism remains 
elusive. 



III. HUBBARD MODEL FOR HYDROGENIC 
SYSTEMS 

A. Overview and formulation 

An important question that can still be asked of a sys- 
tem with positional disorder is whether or not the ground 
state is spin polarized (resulting in macroscopic spin de- 
generacy). In the remainder of this paper, we attempt to 
answer an even more basic question - does there exist, 
even on the nanoscale, large spin degeneracy in systems 
of hydrogenic centers, using an appropriate Hubbard-likc 
description? The paramount conclusion is that there does 
exist a regime in doped semiconductors which is more 
amenable to Nagaoka ferromagnetism. Interestingly, this 
regime is attainable in nanoscale quantum dots and het- 
erostructures, but not accessible to bulk systems. There 
we find Nagaoka-like ferromagnetism in the presence of 
disorder, at least at the nanoscale, and that this regime 
also possesses a higher likelihood of emerging on meso- 
scopic or macroscopic scales (e.g. in modulation doped 
systems). In this section we introduce and motivate 
the generalized Hubbard model used to characterize the 
doped semiconductor problem. 

B. Random Hubbard Model: positional disorder 

As a first approximation, a system of N s randomly po- 
sitioned donors can be modeled with the Hubbard Hamil- 
tonian obtained by adding site-dependence to the hop- 
ping amplitude in Eq. @. Specifically, we make a 
function of the site separation: tij = t(\ri — rj\), result- 
ing in the Hamiltonian: 

Tirdm = - X! (*« C L C Jff + n ' C ') + U n *T n ii ( 9 ) 
i,j,a i 

where i,j = l... N s . 

This takes into proper account the random position- 
ing of the donors, and, as we discuss below, should be 
a good model for both uncompensated and compensated 
bulk semiconductors with < 1 electron per donor site 
(in the latter case, a more rigorous treatment would ad- 
ditionally include random on-site energies reflecting the 
random fields generated by the (positively charged) ac- 
ceptor sites). 

C. Hubbard model generalization: 
occupation-dependent hopping 

A shortcoming of Ti r dm (Eq. ©), both for the lattice 
and random case, is that it does not account for a fun- 
damental property of hydrogen: the two-electron wave- 
function of the H~ ion has much greater extent than the 
one-electron wavefunction of the H atom. This is re- 
flected in the binding energy (the energy required to re- 
move an electron) of H~ being only 0.0555 Ry*, whereas 
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1 Ry* is necessary to remove the electron of H&2^ In- 
deed, using that an effective Bohr radius a* scales as 
1 / ^/-Ebinding, we find that the ratio of Bohr radii for H~ 
and H, a* H -/a* H = \/n)/V0.0555 « 4, showing that the 
wavefunction of H~ is several times larger than that of 
H. Variational treatments of the H~ ion (21 as well as 
an effective pseudopotential calculation,— determine the 
ratio to be in the range 2 — 4. This affects the Hubbard- 
description of the system because it is much easier for 
an electron on a doubly-occupied hydrogenic center to 
hop away than it is for the electron on a singly-occupied 
site to make a similar hop. This implies that the hop- 
ping amplitude seen by an itinerant electron, hopping 
around in a background of singly-occupied sites, is larger 
than that seen by a hole in a similar background. The 
fact that the ratio of the two radii is substantial (2 — 4), 
and the hopping amplitude is exponentially dependent 
on the radius (in the low density regime) , suggest that a 
doped semiconductor above half-filling is in a quite differ- 
ent regime of parameters than the conventional compen- 
sated semiconductor (a system below half-filling) . Such a 
regime, while not obtainable in bulk doped semiconduc- 
tors, should be realizable in semiconductor heterostruc- 
tures, as well as quantum dots. In Hubbard model par- 
lance, near half-filling the hopping amplitude for an elec- 
tron is much larger than for a hole. At the very least, 
the different radii of the doubly- vs. singly occupied sites 
suggest that we modify lattice Hubbard Hamiltonian ((2|) 
to become: 

Ti* = - 53 n 3) c \a c i° + n - c -) + U 53 n ^ Ul i 

(10) 

where rii is the total occupation of site i, and the hopping 
now has occupation dependence given by the piecewise 
function (the hopping corresponding to the different am- 
plitudes t and t is shown pictorially on the right): 

t(rn, rij) = 



t rij — 1 , rii = 2 




t otherwise 




where t is larger (and as we will see, can be much 
larger) than t£^ This model enhances the hopping from 
a doubly-occupied site to an already singly occupied site 
(which will become doubly occupied after the hop) . One 
may question why the hopping from a doubly-occupied 
site to an empty site (the middle picture of Eq. (fTTj) ) 
is not also enhanced. The primary reason is that the 
present formulation is the only way, within the single- 
band Hubbard model, to preserve the asymptotic spa- 
tial dependence of the effective exchange interaction: 



J(r) - e - 2r K (recall J - t 2 /U and t ~ e' r ' a ^). This 
is of essential importance, since this relation for J has 
been shown to be asymptotically exact.— 

Note that Eq. (fT0|) is in general not electron-hole sym- 
metric. Only when t = t and the system is on a bipartite 
lattice is electron-hole symmetry preserved^ The gen- 
eral lack thereof is readily seen, since an itinerant hole 
hops with amplitude t whereas an itinerant electron hops 
with t. Indeed, the effective low energy theory of the 
three-parameter Hubbard Hamiltonian when there is less 
than one electron per site (below half- filling), in the limit 
U t, is independent of t and given by the familiar t — J 
Hamiltonian: 

Htj = -t 53 (i 1 -n i a)4 a Cj (7 (l-n0)+h..c?j 

<ij>cr 

+ J 53 (Si-Sj - -niiij) 

<ij> ^ ' 

where the AF exchange J = 4t 2 /C7, c| CT (ci a ) is the elec- 
tron creation (annihilation) operator, and the spin oper- 
ator Si is as previously defined. When there is greater 
than one electron per site, however, the low energy spec- 
trum (in the large U/t limit) is given by a t — J model, 
where t is replaced by t in Eq. (fT2|) . (l — m a is replaced by 
riin, and where J remains determined by the Hubbard t 
parameter, as one might expect; 95 The Hilbert space re- 
striction then excludes doubly- vacant sites. It is worth 
noting that in the usual t — J model on a non-bipartite 
graph (defined as a set of sites and hopping links), ex- 
cluding doubly-vacant sites is not equivalent to exclud- 
ing doubly-occupied sites. This directly corresponds to 
the lack of electron-hole symmetry in the corresponding 
Hubbard problem. 

It is important to remember that the electron creation 
and annihilation operators in these models act on a sys- 
tem with a fixed number and arrangement of sites. In a 
semiconductor, each site corresponds to a dopant atom, 
and when we speak of adding electrons or holes to the 
system we mean addition or subtraction of carriers while 
leaving the underlying dopant configuration fixed. Thus, 
the electron-hole asymmetry here is not an asymmetry 
between n-type and p-type semiconductors, but an asym- 
metry between a doped semiconductor which has more 
electrons than dopant atoms and one which has less elec- 
trons dopant atoms. 

Hirsch has investigated a similar Hubbard model with 
occupation-dependent hopping, but in a different regime 
with its focus on superconducting pairs^i We proceed 
with semiconductors in mind, and to allow for the ran- 
dom placement of sites, we add positional dependence to 
the hopping amplitude in Eq. (|10[) . similar to the modi- 
fication yielding Eq. ^ earlier, to arrive at: 

(13) 
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where ni is the total occupation of site i, and ijj now has 
an occupation dependence given by: 



One way to view the manifest electron-hole asymme- 
try of models (fT0|) . IT2|) . and (Tl3|) is that systems above 
half-filling are effectively /ess random, and hold greater 
hope for the Nagaoka phenomenon to take place. This 
reasoning follows from electrons having more extended 
wavefunctions than holes, and the concomitant existence 
of two distinct length scales. Because the electron wave- 
functions average over much more of the disorder, sys- 
tems with a small percentage of extra electrons expe- 
rience a greatly reduced effect of the positional disorder 
when compared with corresponding hole-doped (i.e. com- 
pensated) systems, and so behave more like the uniform 
lattice. 

Hope for Nagaoka ferromagnetism in electron-doped 
semiconductors is also found by considering the relation 
of conventional doped semiconductors to diluted mag- 
netic semiconductors (DMS), for which ferromagnetism 
does co-exist with disorder. In one type of DMS (III-V), 
a transition metal atom acts as both a dopant and a local 
moment (coming from the unfilled d-shcll of the atom). 
For instance, in Gai-^Mn^As ) 97 ' 98 the Mn atom acts as 
an acceptor (p-type) and local moment. These systems 
also have substantial disorder (due to dopant positions 
and anti-site defects in the semiconductor itself, e.g. As 
on Ga sites), but possess macroscopic ferromagnetism 
for temperatures up to 100K! 97 Thus, disorder by itself 
does not always destroy ferromagnetism; in fact, in some 
cases it may even enhance the ferromagnetic transition 
temperature i"i 100 

One important difference between conventional "non- 
magnetic" doped semiconductors and DMS is that there 
exists in the latter two different length scales - the Bohr 
radius of the Mn hole wavefunction (~ 10 A) and the 
extent of the localized spin on the Mn (~ 1 — 2 A) . Thus, 
each hole's wavefunction extends over several Mn spins, a 
phenomenon which is only accentuated as holes delocal- 
ize further at higher Mn density. This allows the carrier- 
magnetic moment interaction to dominate, resulting in 
a FM ground stat o 99 ' 101 . In the electron-doped semi- 
conductor, the Bohr radius of the electrons that singly- 
occupy sites (which give rise to the effective AF exchange 
interaction J ~ t 2 /U), is much smaller than the radius 
of the electrons which doubly-occupy a site. This di- 
chotomy of length scales could similarly conspire to re- 
sult in carrier hopping being dominant and ultimately a 
ferromagnetic (Nagaoka) ground state. (The other dif- 
ference, of course, is the existence of multiple bands in 
DMS, which facilitates FM.) 



D. Parameter ranges and calculation details 

The first step in our analysis of the Hamiltonians (jTUJ) 
and (TT3|) is to find values (or ranges of values) appropri- 
ate for their parameters. The models are described by 
the dimensionless ratios t/t and U/t (which depend on a 
pair of site indices in the case of Eq. (|13p). To find values 
of U ft and t/t appropriate for doped semiconductors, we 
performed a calculation of the single particle states of 
donors placed on a simple cubic lattice. Note that al- 
though much of our work deals with 2D systems, atomic 
hydrogen is intrinsically a 3D problem, and thus the cal- 
culation of realistic parameters for a system of many hy- 
drogenic centers should likewise be in three dimensions. 
We choose the simplest such 3D arrangement of centers, 
the simple cubic lattice. 

As already stated, a hydrogen (H) atom binds its elec- 
tron with a strength of 1 Ry* and will bind a second 
electron with 0.0555 Ry* to form a H~ ion. If all of the 
dopants are positioned on a superlattice, then these two 
levels broaden in the usual manner into two impurity 
bands. The exact details of the particle bands depend 
on the spin configuration in the ground state. Due to 
the H~ ion's wavefunction being more spatially extended 
than that of the H atom, the width of the upper impurity 
band is significantly greater than that of the lower band. 

We have calculated these bands for a ferromagnetic 
configuration of spins in the ground state of a filled lower 
band (i.e. the uncompensated case). We follow Bhatt 
and Rice^ and use pseudopotentials and a sphericalized 
Wigner-Seitz (WS) method on a cubic superlattice. De- 
tails of the band calculation can be found elsewhere.— We 
then extract the dependence of t and t on the impurity 
density (or equivalently, on the lattice constant) by fit- 
ting the calculated bandwidths to a tight binding model. 
Using the well-known tight binding relationship between 
hopping parameter and bandwidth on un-frustrated lat- 
tices yields (where z is the lattice coordination number): 

2zt = width of lower band 
2zt = width of upper band (15) 
U = band gap at zero density . 

We find U « lRy* and, by matching the bandwidths for 
the 3D case, we obtain the tight binding parameters t(b), 
t(b). Figure [5] shows the dependence of the dimension- 
less Hubbard parameter ratios on the superlattice spacing 
(lower axis) and impurity density (upper axis). It shows 
clearly that the range of U/t and t/t can be varied sub- 
stantially in the doped semiconductors. The large span of 
U/t originates in the exponential dependence of the hop- 
ping parameter on the atomic spacing, and the variation 
of t/t from the relatively large size of the two-electron 
wavefunction appearing as a factor in this exponential. 

In the results that follow, we either use the exact pa- 
rameter ratios found here or consider the effect of vary- 
ing the parameter ratios within the ranges U/t = [5, 100] 
and t/t = [1, 10], which are conservative when compared 
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ing (related to the dopant density p by p = -jkg , so the metal- 
insulator occurs at R c /aB = 4). 



ditions. Note that only nearest neighbor links are kept 
in the model (see Eq. (fTU)) ), so that there is a single pair 
(t, t) of kinetic parameters. We refer to a lattice as being 
bipartite or non-bipartite if the corresponding Hubbard 
model with only nearest neighbor hopping is respectively 
bipartite or not. Section [IV Bl presents results from clus- 
ters with open boundary conditions and selected struc- 
tures for which all nearest neighbors are equidistant (so 
there is again a single pair of kinetic parameters) . We use 
the term cluster in this section to refer to a finite system 
possessing less symmetry than a finite lattice. In section 
IIV C[ clusters constructed to have only two or three pairs 
of kinetic parameters are considered with open boundary 
conditions. There we also describe a method of adding 
random perturbations to clusters, and present results for 
several cases. Finally, in section TlV Dl we analyze ensem- 
bles of random clusters. We generate these ensembles 
with a fixed density, and exact diagonalization results of 
the individual clusters are averaged to produce our final 
results. 



to the physically attainable ranges. After determining 
the parameter ranges of interest, we solve both Hubbard 
and t — J models on finite systems. We numerically find 
the ground state, and determine how its spin depends 
on t/t, U/t, system size, and system geometry. Hamil- 
tonians (flQ|) and (fl~3|) were solved using exact diagonal- 
ization, ultimately using a generalization of the Lanc- 
zos method*^ With four states allowed on a site, the 
Hilbert space grows exponentially in the number of sites, 
restricting the size of tractable systems significantly. Sev- 
eral optimizations have been exploited to push back this 
computational barrier. First, since both Hubbard and 
t — J Hamiltonians commute with the z-component of 
total spin, it follows from the properties of the SU(2) 
group, that we can restrict the the Hilbert space to the 
minimal S z sector without reducing the support of the 
spectrum. Second, all spatial symmetries are utilized via 
group theoretic techniques to divide the Hilbert space 
into sectors for which the Hamiltonian matrix is block 
diagonal. Third, we factorize the action of the Hubbard 
Hamiltonian into "up spin" and "down spin" parts, al- 
lowing more efficient computation of the matrix elements. 
In the t — J model, this can be done only for the kinetic 
term. 



A. Finite Lattices 

We have solved the nearest-neighbor Hubbard and cor- 
responding t — J models on finite square (8, 10, and 16 
sites), honeycomb (6 and 10 sites), and triangular (7 and 
9 sites) lattices. These are shown in Fig. [5] with the sites 
of a single unit cell connected, so that the method of ap- 
plying periodic boundary conditions in each case is clear. 
Note the choice of unit cell for all of the bipartite lattices 
(square and honeycomb) allows a classical Neel state spin 
assignment, where all of a site's nearest neighbors have 
spin opposite to it. This requirement is important since a 
finite bipartite lattice that is magnetically frustrated due 
to boundary conditions may have an exaggerated prefer- 
ence for FM. 

Each finite lattice, with periodic boundary conditions, 
was doped with up to two electrons or holes away from 
half-filling. Denoting the number of electrons N e , this 
means that N s - 2 < N e < N s + 2. The Hubbard 
model depends on the two dimensionless ratios U/t and 
t/t, whereas the t — J model depends only on t/J = 
j(t/t)(U/t). Thus, the value of t/J marking the onset 
of the Nagaoka state defines a straight line in log U/t 
vs. logt/t space with slope —1. We consider each lattice 
in turn below. 



IV. RESULTS FOR GROUND STATE SPIN IN 
FINITE CLUSTERS 

Here we present the results of solving our generalized 
Hubbard model on finite systems. The results and dis- 
cussion are divided into units based on the amount of 
structure present in the system, and what type of bound- 
ary conditions were used. Section TlV Al considers systems 
with finite lattice structure and periodic boundary con- 



1. Square Lattice 

The square lattice, the stereotypical 2D lattice, is bi- 
partite and is itself a Bravais lattice. Figure [7] shows 
the ground state spin phase diagram for the 8-, 10-, and 
16-site square lattices doped with one electron, up to 
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FIG. 6: Lattice geometries for the square, honeycomb and triangular lattices used in this section. The lines connect sites of 
the finite lattice, which is repeated to show how periodic boundary conditions are implemented. 



t/t = 5. One sees that an increase in t/t causes the re- 
gion where the ground state attains its maximum spin to 
increase. This confirms our intuition about the model, 
that a FM ground state is more likely when the carri- 
ers (an extra electron in this case) have greater hopping 
amplitude. (Recall that a greater hopping amplitude in- 
creases the kinetic energy gain of a delocalized electron 
in a background of aligned spins relative to the case when 
the background spins are in an AF or random arrange- 
ment.) Up to t/t = 5, the minimal U/t needed for a fully 
polarized ground state falls roughly as a power law with 
t/t. The i — J model gives a fairly accurate fit to the 
Hubbard data (predicting a power law with exponent -1, 
shown by the lines in Fig. [7]). The fit is especially good at 
low t/t, which coincides with larger U/t values and thus 
is where we expect the t— J model to be most accurate. 
Beyond t/t = 5, the same general trend is observed, but 
the phase diagram becomes more complicated as regions 
of intermediate polarization arise, making the transition 



from low spin to maximal spin less abrupt (and closer 
to a second order transition). This behavior is shown in 
Fig. [5] for the 16-site square lattice with N e = 17. The 
t — J line in this case runs through the regions of in- 
termediate spin, though the model itself gives a direct 
transition from minimal to fully saturated ground state 
spin. 

A comparison of these electron-doped systems with 
corresponding hole-doped systems reveals a pronounced 
electron-hole asymmetry. This is expected from the 
model, since for t ^ t the Hamiltonian is not electron- 
hole symmetric: electrons hop with t whereas holes hop 
with amplitude t. Figure [5] compares the Hubbard model 
with N e — N s ± 1 (one extra electron or one hole) on fi- 
nite square lattices. In the larger 10- and 16-site lattices 
with one hole we see very little dependence of the ground 
state spin on t/t, as would be naively expected. [In the 
8-site square lattice an increase in t/t actually hinders 
ferromagnetism, seen by an increase in the U/t neces- 
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FIG. 7: (Color online) Ground state spin diagram resulting 
from the exact diagonalization of Eq. (|10[) on 8-, 10-, and 16- 
site square lattices (periodic b.c.) with 9, 11, and 17 electrons 
respectively. Hubbard model results are displayed as open 
symbols. Lines show the result of the corresponding t — J 
model as described in the text. S max denotes the region of 
largest allowed spin (actual value depends on the lattice size), 
and Si ow marks the region of unsaturated (usually minimal) 
ground state spin. 
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FIG. 9: (Color online) Ground state spin diagram for the 
8-, 10-, and 16-site square lattices showing the asymmetry 
between doping with a single hole (dashed line) and a single 
electron (solid line). 




t/t 

FIG. 8: Detailed ground state spin diagram of the 16-site 
square lattice with 17 electrons. Labels indicate the ground 
state's total spin. As t/t increases beyond 3, the transition 
to a maximally polarized state is less abrupt and regions of 
partial spin polarization exist. We find that the t — J model 
gives a direct transition from S = | to 5 = -y-, which is 
shown as the dashed line. 



sary to reach the totally spin-polarized state. This is 
most likely a finite size effect, but may have interesting 
ramifications in the context of finite clusters (see section 
HVCI below)]. It is clear that the asymmetry between the 
electron- and hole-doped results originates from the elec- 
tronic states having greater radius than the hole states, 
since for equal radii (t — t) the square lattice is bipar- 
tite and the problem is electron-hole symmetric. Figure 
[9] is the first of many that illustrate a central result of 
this thesis: high-spin ground states are attained at much 
lower U/t in the electron-doped case than in the hole- 
doped case. 



The ground state spin of the Hubbard model on finite 
2D square lattices with periodic boundary conditions is 
knowr>i2^ to behave somewhat erratically as a function 
of the number of electrons (N e ), and techniques involv- 
ing an average over varied boundary conditions have had 
some success as smoothing out, as well as explaining this 
behavior fi2^ We do not address these issues here; instead 
we focus on the square lattice at two dopings that are 
known to give high-spin ground states when used with 
periodic boundary conditions. In addition to the single 
electron or hole configurations already described, the 16- 
site square lattice with 4 electrons (N e = 20) is known 
to have a ground state spin of maximal value (S = 5). 
Figure [POl shows the effect of varying t/t in this case, and 
we see, similarly to the case of a single carrier, that in- 
creasing t/t decreases the value of U/t needed to attain 
the fully saturated ground state. 
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FIG. 10: Ground state phase diagram of the 16-site square 
lattice with 4 electrons above half-filling (20 electrons total). 
The line is a spline fit, and is provided as a guide for the eye. 



2. Honeycomb Lattice 
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FIG. 11: (Color online) Ground state spin diagram from the 
exact diagonalization of 6- and 10-site honeycomb lattices 
doped with a single electron (i.e. with 7 and 11 electrons 
respectively) showing the boundary of the region where there 
is a complete spin polarization. In the 6-site lattice the tran- 
sition is from S=5/2 to S=3/2, whereas in the 10-site lattice 
the transition is more abrupt, changing from S=9/2 to S=l/2 
within the resolution used. 



There has been a revived interest in the honeycomb lat- 
tice since the recent surge in graphene-related research. 
Though it is not itself a Bravais lattice (it is a triangu- 
lar lattice with a two-point basis) , the honeycomb lattice 
is bipartite and thus the Hubbard model is electron-hole 
symmetric on it for t — t. The mean-field ground state 
phase diagram of the t = t Hubbard model for hole- 
doped systems shows the existence and stability of the 
Nagaoka phase at large U/t near half-fillingji2^ The mag- 
netic ground state diagrams for Hamiltonian (|10p on 6- 
and 10-site honeycomb lattices with one electron or hole 
away from half-filling (N e — 2V a ±l) are shown in Figs. [Til 
and [12] respectively. 

We find similar qualitative behavior to that of the 
square lattices: for systems with N e = N s + 1, increasing 
t/t expands the region of phase space for which the spin 
is maximal. Again, the t — J model result agrees well 
with the Hubbard results for low t/t. In the case of sin- 
gle hole-doping (N e = N a — 1), there is little dependence 
on t/t in the 10-site lattice whereas there is the opposite 
t/t dependence in the smaller 6-site lattice, similar to the 
case of the 8-site square lattice. 



3. Triangular Lattice 

The triangular lattice is a Bravais lattice of particular 
interest, since it magnetically frustrated (not bipartite). 
A recent study of the triangular latticed using a many- 
body expansion technique finds that, at large U/t, a 
120 °-ordered AF phase is stable at and below half-filling, 
and becomes unstable above half-filling. In past studies 
of finite clusters, it was likewise found that at half-filling 
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FIG. 12: (Color online) Exact diagonalization results showing 
the boundary of the fully spin polarized region on the 6- and 
10-site honeycomb lattices doped with a single hole (i.e. with 
5 and 9 electrons respectively). In the 10-site case, the spin 
on the unsaturated side of the transition is S — | except for 
a region of S = | found at intermediate U/t for t/t > 10; 
on the 6-site lattice the unsaturated state has uniform spin 
|. Note that there is much less variation with respect to t/t 
when compared with Fig. 1111 
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antiferromagnetic states are optimal in non-bipartite sys- 
tems (due to the quantum fluctuations arising from what 
would be frustrated bonds in a static picture).— 

With a single extra electron (N e = N s + 1), the Hub- 
bard model on 7- and 9-site lattices displays saturated 
ferromagnetism very strongly (on the 9-site lattice with 
t = t, U/t « 15 results in a spin polarized ground state). 
Figure [13] shows our results for the Hubbard model on 
finite triangular lattices with one extra electron. Classi- 
cally, the observed dominance of ferromagnetism could be 
linked to a suppression of competing AF configurations 
(frustrated on the triangular lattice). One must be care- 
ful, however, when applying this reasoning to quantum 
models, as studies have shown that antiferromagnetism 
is enhanced on the triangular lattice with a single hole- 2 
due to the subtle interplay of quantum phases. The reg- 
nancy of ferromagnetism may also be due to the large 
number of tight loops in the lattice. Pastor et aZ.— have 
remarked that the presence of triangular or square loops 
coincides with ferromagnetism in finite clusters, and we 
reach similar findings in our study of clusters below (see 
sections IIVBI and IIV C|) . The strong FM we see here 
suggests that this connection extends to lattices as well. 

The t — J data for the triangular lattice fits the Hub- 
bard data less well than in the previous bipartite lattices. 
For the 9-site triangular lattice the t — J result underes- 
timates the region of saturated spin, and in the case of 
the 7-site triangular lattice, the Hubbard model does not 
even transition to the unsaturated state predicted by the 
t — J model. The discrepancy is not an immediate cause 
for concern, and might even be expected, given the low 
U/t values at which the the transitions occur. 

Since the triangular lattice problem is not bipartite, 
there can be (and is) electron-hole asymmetry even when 
t = t. Figure [14] shows the ground state phase dia- 
gram for single hole-doped 7- and 9-site triangular lat- 
tices (N e = N s — 1). These plots are qualitatively dif- 
ferent from those of the the hole-doped square and hon- 
eycomb lattices: the high-spin region is unsaturated and 
lies at lower U ft than a minimal-spin region which dom- 
inates at large U/t. As t/t is increased, the partially 
polarized region expands up to larger U/t values. The 
mechanism for this may be related to the "kinetic antifer- 
romagnetism" studied by Haerter and Shastry, 72 which 
explains how the phase dependence of a single hole's mo- 
tion enhances antiferromagnetism. 



B. Selected Symmetric Clusters 

Next we consider a select group of two-dimensional 
Hubbard clusters that, like the finite lattices, have only 
a single pair of hopping amplitudes, t and t. Unlike the 
lattices, these clusters are given open boundary condi- 
tions. This corresponds to the physical situation in which 
a small number of sites (dopants or quantum dots) are po- 
sitioned in a plane such that every pair of nearest neigh- 
bors is equidistant. Pastor et a^ 61 > 63 have studied the or- 
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FIG. 13: (Color online) Ground state spin diagram from the 
exact diagonalization of 7- and 9-site triangular lattices when 
doped with a single electron, showing the region of saturated 
spin. On the 9-site lattice, the unsaturated region is pre- 
dominantly 5 — except for a sliver of 5 — 2 close to the 
transition. There is no transition on the Hubbard 7-site lat- 
tice, which has a maximally polarized ground state (5 = 3) 
for the entire plotted area. In the corresponding t ~ J model, 
however, the 7-site lattice has a transition from 5 = 3 to 
5 = 2 near t/J m 3.0 (shown by the dotted line). 
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FIG. 14: (Color online) Ground state spin diagram for the 
7- and 9-site triangular lattices doped with a single hole. 
Nowhere is the ground state spin saturated. Instead, there 
is a region of minimal spin (5 = 0) at large U/t which is en- 
croached upon by a region of partial spin polarization (5 = 2 
and 5 = 3 for 7- and 9-sites respectively) as t/t increases. 
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dinary Hubbard model (Eq. {2])) on all possible geomet- 
rically realizable clusters in two and three dimensions. 
Our analysis of cluster structure here is not as exhaus- 
tive, but we calculate the phase diagram along the t/t 
axis. Clusters are chosen to lie in the plane such as to 
retain some spatial symmetries, and their ground state 
spin is calculated for 1 < t/t < 10 and 5 < U/t < 100 
when doped with 1 or 2 electrons away from half-filling 
(in either direction). Figure [T5l summarizes the results, 
giving each cluster's geometric structure and its maxi- 
mal spin as a function of doping. We see that in most 
cases, the highest spin is attained when doped with a sin- 
gle electron, following our expectation that a low density 
of extra electrons will favor spin polarization. Indeed, 
clusters 1-4, 6, and 7, attain their maximal ground state 
spin when doped with one electron. In contrast, clusters 
5 and 9 have greater spin polarization below half-filling, 
and that their polarization is maximal when doped with 
two holes. 

Clusters 1-3, 5, and 7 we call "ring-like", since each 
is equivalent to a 1-dimensional chain of sites with peri- 
odic boundary conditions. In the pair and triangle (clus- 
ters 1 and 2), the spin listed in Fig. [T5]is the only spin 
found in the considered parameter range. Figure com- 
pares the ground state phase diagrams of the remain- 
ing clusters with \N* — N s \ electrons above and below 
half-filling, for all values of N* such that the resulting 
phase diagrams are non-trivial (have at least two spin 
regions). We see from Figs. [K] and [TH] that the trian- 
gle and square show the greatest percentage spin polar- 
ization above half-filling (both have maximally polarized 
ground states; for the square at large t/t). Also note 
the t/t dependence of the square with one hole vs. with 
one electron, where we see behavior similar to that of 
the square lattices. The pentagon is unusual in that it 
has higher ground state spin when hole-doped. A fully- 
polarized ground state occurs at large U/t when the sys- 
tem is doped with two holes. Lastly, the hexagon shows 
very little t/t dependence, though with two holes (4e~) 
larger t/t creates an interval in U /t with low spin (S — 0). 
This behavior was also seen in the hole-doped bipartite 
lattices of section ITV Al 

The remaining (non-ring-like) clusters, 4, 6, 8, and 9 of 
Fig. [m are created by adjoining triangles and squares. 
This was done with the hope of engineering clusters with 
a high-spin ground states, given the individual properties 
of the triangle and square. Detailed ground state phase 
diagrams for these clusters are presented in Appendix [Al 
We see in general that increasing t/t enlarges the high- 
spin region of the phase diagram for electron-doped clus- 
ters and, in this sense, indicates that the high-spin state 
has become more robust. In hole-doped systems we see 
a much weaker dependence on t/t, and in the clusters 
8 and 9 we see the opposite behavior: as U/t increases 
there is a transition to lower ground state spin. Upon 
electron-doping, we find a correlation between structures 
that have a large number of triangular or square loops 
and those with high spin ground states. This relationship 




FIG. 15: Summary of clusters that have a single pair of 
hopping parameters. S^ M is the maximum spin obtained in 
the window t/t £ [1, 10], U/t S [5, 100] when the system has 1 
or 2 holes or electrons away from half-hlling (x = lh, 2h, le, 2e 
respectively). Note the correspondence of high-spin states 
with larger numbers of tight loops. 



has also been seen in previous worki^ Though a precise 
reason for this correspondence has not been found, we 
believe it is due to such systems being electronically un- 
frustrated, allowing an electron to easily hop among all 
the sites and to be very effective at increasing the ki- 
netic energy of the FM state. Whatever the mechanism, 
a heuristic rule for constructing clusters with high spin 
ground states is that a large cluster with many tight loops 
(triangular or square) is likely to be strongly magnetic. 
This has recently become relevant to experiment through 
the work of Schofield et aL^S. who are able to position 
phosphorous dopants within bulk silicon to nanometer 
accuracy using a scanning tunneling microscopy (STM) 



15 



Geometry 




Ground state phase diagrams 



25 
20 

7 10 



s = 1.1 



3e" 



V 



S = 0.5 



1.5 2 2.5 3 

t/t 




25 



15 



10 - 









3e" 




s = 


1.5 














s= 


0.5 





10 



t/t 



5 = 0.5 



3 5 7 

t/t 



10 



10 - 



5=1 



S = 0- 



U 



10 



8e" 



10 



t/t 



S= 1 



t/i 



10 



FIG. 16: Ground state spin (T = 0) phase diagrams in the U/t — t/t plane for clusters 3, 5, and 7 from Fig. 1151 These 2D 
clusters are "ring-like" in the sense that they are equivalent to ID chains with periodic boundary conditions. The fixed electron 
number is given in the upper-right corner of each plot, and only selected non-trivial diagrams are shown. 



tip. Such capability allows for the construction of cluster 
geometries made "to order" , and opens an entirely new 
area of application for our work. In particular, the ability 
to test for FM behavior (i.e. high spin ground states) in 
finite lattices of dopants would be very valuable. 



C. Distorted clusters 

More complex 2D clusters are obtained by allowing 
more than one pair of hopping parameters (i.e. hop- 
ping is allowed between sites of different separation 
distances). In this section we consider clusters with 
two and three pairs of distinct hopping parameters 
Uti,ti) : i G (1,2,3)}. Some of these can be viewed as 
geometric perturbations of clusters in the last section, 
while many are new geometries not possible under the 
restriction of equidistant nearest neighbors. For a select 
group of clusters with two pairs of hopping parameters, 
we consider the ground state spin as a function of t2/t\ 
and U/ti at a uniform fixed ti/ti, i — 1,2. Our anal- 



ysis is done over the substantial region of phase space: 
h/h S [1, 10], h/U S [0.01,0.5]. (Note that this extends 
to U/t < 10, outside the physical range found earlier, but 
in the direction that favors non-ferromagnetic behavior.) 
The results are summarized in Fig. 1171 which show for 
each geometry the maximal spin achieved with a doping 
of up to two electrons or holes (the maximum is taken 
over the region of phase space stated above). Again we 
find that most clusters attain their highest spin when 
doped with le" (clusters 1, 2, 4, 7, 10, 12, 14, 15, 18, 20, 
and 22). Some of the larger clusters also have high spins 
when doped with two electrons (clusters 11, 18, 20, and 
23), since their density is still low enough to favor FM. 
Although in most cases the maximal spin is greater for 
electron-doping than hole-doping, there are some which 
attain high spins even when hole-doped (e.g. clusters 8, 
9, 11, and 15). 

We focus on the ground state spin behavior of three 
clusters from Figs. [TTJ 11, 12, and 20. Ground state 
phase diagrams showing the spin for these clusters are in 
the Fig. [TS1 Each row of the table shows the geometry 
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FIG. 17: Summary of maximum ground state spins for clusters that have two pairs of kinetic parameters (two distinct nearest 
neighbor distances). Solid lines represent hopping amplitude ti, and dashed lines t2- Cluster geometries are listed by size, and 
maximal spin is given for dopings of -2,-1,1, and 2 electrons away from half-filling. Each cluster is identified by a number, # c ;, 
and the maximum is taken over the region t-2.jt\ € [1, 10], ti/U G [0.01, 0.5] for ti/U uniformly set = 1, 5, and 10. 



and two ground state phase diagrams of a cluster with 
a fixed number of sites N s and electrons N e . The two 
diagrams correspond to t/t = 1 and 5, as indicated by the 
column headings. The charge of the cluster Q — N s — N e 
(the negative of its doping relative to half filling) is given 
in the third column. For each selected cluster, phase 
diagrams are only shown for Q = ±1. The transition lines 
in these plots are found by finding the ground state spin 



on a grid in parameter space, then fitting the transitions 
between grid points with smooth curves. Detailed phase 
diagrams of all non-trivial cases are given Appendix [B] 

Clusters 11, 12, and 20 have fully spin-polarized 
ground states when doped with one electron, and for 
this reason will be used as starting points in later per- 
turbation schemes. The movement of ground state spin 
boundaries as t/t is increased in steps (t/t = 1,5, 10) is 




FIG. 18: Ground state spin diagrams for selected clusters from Fig. 1171 
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seen in each row of the table. In cluster 11 the region of 
t 2 /ti vs. U jt\ space with maximal spin expands for both 
electron- and hole-doped cases as t/t increases, which is 
interesting since the effect of a larger t/t on a hole-doped 
system is expected to be relatively minor. In cluster 12, 
a similar increase in polarization with larger t/t is only 
seen in the single electron-doped case (the Q = — 1 case 
is all that is shown, since all other dopings have minimal 
spin throughout the plotted region; see Fig. [17]). Clus- 
ter 20 behaves very much like we naively expect: in the 
hole-doped case (Q = +1) the diagram is almost insen- 
sitive to changing t/t, while for Q = —1, the region of 
maximal polarization clearly expands at the expense of 
other lower-spin regions. We note that in all cases high- 
spin ground states occur when t 2 /ti is close to 1, that is, 
when the dotted hopping links in the tables are nearly as 
strong as the solid links and the triangles and pairs that 
make up cluster are more strongly coupled. 



Several major conclusions may be drawn from this 
data. First, there are many instances of high-spin 
ground states among these clusters, many of which can 
be thought of as a weak coupling fa) between triangles 
and pairs with a stronger internal coupling fa). In a 
real system, where the broad distribution of inter-site 
distances due to positional randomness creates exponen- 
tially strong and weak bonds, these results give some 
hope that the spin-polarization seen in the isolated tri- 
angle, for example, will survive in the presence of pertur- 
bation due to other sites, and that this interaction may 
even lead to spin polarization on longer length scales. 
Second, it is found almost universally that increasing t/t 
leads to greater spin polarization in electron- doped clus- 
ters, just as in the finite lattices ( section HV A[) and single- 
hopping parameter clusters (section IIVB[) . In electron- 
doped clusters, we continue to see a correlation between 
the number of triangular loops in a cluster and that clus- 
ter's maximal spin. For instance, compare clusters 5 and 
7 with clusters 14 and 15 of Fig. [T7| (the latter are much 
more magnetic) . In hole-doped systems we generally find 
lower spin values, and often there is a high-spin region 
at low U jt\ . This inverted relationship in clusters below 
half- filling was also found in section IIVBI and on the 8- 
site square lattice. Lastly, we note that although there 
is potential for high-spin states, there are many clusters 
that have large regions of minimal ground state spin. We 
find overall that the Nagaoka-like ferromagnetic effect we 
observe is very sensitive to geometry, though the sensi- 
tivity decreases at large t/t. 



Next, we test the stability of a select few of the high- 
spin ground states found above. For clusters 11, 12, and 
20 of Fig. Q~7] we further reduce the spatial symmetry 
by additional geometric distortion, as shown in Fig. [T9j 
The distortion introduces a third pair of hopping ampli- 
tudes fa,ts), and the ratio £3/^2 measures the amount 
of distortion. 
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FIG. 19: Geometries of clusters obtained by geometric distor- 
tion of clusters 11, 12, and 20 of Fig. [17] with three pairs of 
kinetic parameters t\ > ti > £3. 
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FIG. 20: Ground state spin diagram for distorted clusters. 



We fix tz/ti at a value for which the undistorted 
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(*3 = h) cluster has a high-spin ground state, and de- 
termine the amount of distortion that can be applied 
(i.e. the lowest value t^jt^ can attain) before the clus- 
ter loses its high spin state. The value of ti/U is fixed 
(i.e. in each run, all of the links forming the cluster have 
the same t/t ratio), and the resulting ground state phase 
diagrams as a function of £3/^2 and U/t\ are shown in 
Fig. 1201 There are two key points resulting from this 
data. First, as t/t becomes larger, the high-spin ground 



states become more robust to the geometric fluctuation 
considered here: regions with high-spin ground states 
persist to lower values of t^/t% as t; L jt; L is raised. (Recall 
that lower t^/t^ corresponds to larger geometric distor- 
tion.) Second, the high-spin ground states are more ro- 
bust at larger U/ti, since the curves for fixed U/ti move 
to lower values of t%/t\ as U increases (e.g. U — 100 curve 
lies below the U = 50 and U = 20 curves). 
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FIG. 21: (Color online) Result of randomizing clusters 11, 12, and 20 of Fig. fT71 For clusters 11 and 12, U/ti = 20 and 
£2/^1 = 0.1; for cluster 20, U/ti = 100 and t%/t\ = 0.3. We set U/ti to 1, 2.5, and 5 as indicated by the row headers. 



To further probe the robustness of a given cluster's 
high spin ground state, we consider multiplying each tj 
of the cluster by a random factor A whose logarithm is 
chosen from the box distribution P(logA) = l/(21oga), 
log A € [— log a, + log a) where a > 1. Thus, when a = 1 



the system is unperturbed, and for a > 1 each hop- 
ping amplitude is independently multiplied by a different 
random number between 1/a and a. Compared to the 
specific geometrical distortions analyzed in the preceding 
paragraph, this method of introducing randomness more 
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accurately characterizes the fluctuations we expect in a 
real system, since the hopping is exponentially dependent 
on the inter-site distance and we do not expect the fluc- 
tuations to preserve any symmetry present in the cluster. 
We take as our starting point a cluster known to have a 
high spin ground state and average over 1-5 thousand of 
the just described random perturbations. Then, we tab- 
ulate the percentage of the randomly perturbed clusters 
possessing each possible value of the ground state spin. 
The shaded regions in plots of Fig. [5T] show how these 
percentages vary as a function of a, with the different 
figures corresponding to initial clusters 11, 12, and 20 of 
Fig. [lTj The boundaries of the regions are spline fits to 
the data. We set U/t — 20, a relatively low value for 
doped semiconductors, to more clearly see the effect of 
a (at larger U/t the high-spin ground state becomes in- 
creasingly robust). We see in all cases the general move- 
ment, in a probabilistic sense, of the clusters from high 
to low spin as a is increased, but that this effect is signif- 
icantly mitigated by raising t/t. As t/t becomes larger, 
the percentage of the clusters that retain the high-spin 
ground state of the original (a = 1) cluster grows sub- 
stantially. Thus, we again find that increasing t/t makes 
high spin ground states significantly more robust to ran- 
dom geometric fluctuations, this time to fluctuations sim- 
ilar to those we expect in an actual doped semiconductor. 
This result gives additional hope for the viability of con- 
structing magnetic clusters using an STM tip (described 
above), where there would inevitably be slight errors in 
the dopant positions. 

D. Randomly distributed finite clusters of fixed 
density 

In sections IIVBI and IIV CI we solved generalized Hub- 
bard and t — J models on a variety of clusters that 
were constructed to have some spatial symmetries and 
at most a few pairs of hopping parameters (ti,U). This 
section and the next give an extensive analysis of clus- 
ters with completely random structures and several types 
of boundary conditions. Also, instead of considering a 
range of t/t values, we use only the parameters given by 
our realistic band calculation described in section UlIDI 
In (^-dimensions, clusters with N s sites and fixed den- 
sity p are generated by randomly placing N s sites within 
a d-dimensional hypercube of side length L such that 
p(a^ )~ d = N s /L d . We fix U = 1 Ry* and determine the 
hopping parameters fy by setting fy = £ ( | — fj-|), where 
t(r) is given by the lattice calculation described earlier 
(see Fig. [5]). We consider three different models, each 
corresponding to a different method of setting ty : 

1 ■ tij — tij • 

2. Analogous to , using t(r): tij = t(\fi — fj\), where 
t(r) is obtained from the broadening of the upper 
impurity band, referred to as the D~ band in semi- 
conductor literature. 



3. Set t^ = C, where C is a constant. The value of C 
is chosen to be U/2. 

The first case is the regular Hubbard model for ran- 
domly distributed sites, and does not take into account 
the special property of hydrogenic centers. Model 2 
takes into account the larger extent of the D~ state. 
Model 3 is to simulate a situation when the radius of 
the D~ state becomes very large, to see how big an ef- 
fect that would have on the possibility of Nagaoka ferro- 
magnetism. We choose C — U/2 since this is close i(r) 
when r = Og, the smallest separation for which the tight 
binding model could apply. Since t(r) increases with de- 
creasing r, t(a^) w U/2 is of order the maximal t found 
in the entire system. 

Given a fixed cluster size and density, we exactly solve 
many (between 10 4 and 10 6 ) clusters and construct a 
histogram of ground state spin values. Results obtained 
using each of the three models are compared to assess 
the effect of the nature of the doubly-occupied state. We 
have calculated the spin histograms for clusters in two 
and three dimensions with sizes from N s = 4 — 7 and 
for densities p = j^qq, ygg, and in 2D, (correspond- 
ing to ~ 0.005, 0.05, and 0.15 times the Mott metal- 
insulator transition density) and p = g^g, g^, and g|jj 
in 3D (corresponding to 0.01, 0.1 and 0.3 times the Mott 
density). Further, we have considered open as well as 
periodic boundary conditions. In an actual macroscopic 
sample, clusters will be connected to other clusters of dif- 
ferent local densities. Thus, the physical situation will be 
intermediate between the cases of open (where each clus- 
ter is surrounded by no others) and periodic (where each 
cluster is effectively surrounded by others of the same 
density) boundary conditions. The latter is closer to the 
actual case at high densities, the former at low density. 

1. Hopping set by band calculation: tij — t{rij) 

Here we present results for two- and three-dimensional 
random clusters. Clusters have all inter-site links present 
(i.e. hopping is not restricted to be between nearest 
neighbor sites only). We find the distribution of ground 
state spin values for ensembles of clusters with fixed size 
N s , density p, doping (either one extra electron or one 
hole) , and model for determining t^ . 

Raw spin distribution data, shown by tables contain- 
ing the percentage of clusters with each possible spin, are 
given for two-dimensional clusters, with open and peri- 
odic boundary conditions, in Appendix [Cl Correspond- 
ing results for three-dimensional clusters could not be 
included in this pap er du e to length considerations, and 
can be found in Ref. Il07l 

Here, we summarize the data by plotting the average 
spin and the percentage of magnetic clusters (those with 
greater than minimal spin) as a function of doping (zero 
doping = half-filled). We show only the results for 2D 
clusters with open boundary conditions; similar plots for 
periodic boundary conditions can be found in Appendix 
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[Cl Figure [22] shows the average spin of such clusters. 
There is some variation in the average spin due to even- 
odd asymmetry: clusters with an odd number of elec- 
trons have minimum spin S m i n — \, while those with 
an even number have S m i n = 0. To remove this effect, 
Fig. 1231 shows the average spin relative to S m i n (*-e- 0.5 
is subtracted from cases of odd electron number) . A sec- 
ond measure of a systems magnetic behavior is the per- 
centage of clusters with above minimal spin. We define 
any cluster with greater than minimal ground state spin 
(equivalently, spin > 1 since the minimal spin is either 
or 1/2) as a magnetic cluster, and Fig. 1241 shows this 
quantity as a function of doping for the different cluster 
sizes (2D clusters with open b.c). Although both the 
average spin and percentage of magnetic clusters provide 
less detailed information than the spin distribution data 
( Appendix [Cj) . they also suffer less from finite size effects 
and give a more concise picture of the results. 

The study of Figs. [22ll24l and the more extensive data 
of Appendix [Cl and Ref. Il07l reveals several trends. First, 
clusters with periodic boundary conditions tend to have 
a larger total spin than those of the same size and den- 
sity but with open boundary conditions (see Appendix 
[Cl Table IVI[) . This is especially true for large clus- 
ters [N s — 6, 7) and at lower density. This may be 
due to the increased connectedness of clusters with peri- 
odic boundary conditions compared to those with open 
boundary conditions. In a system that is more connected 
(i.e. where there are more nonzero hopping amplitudes 
hj), electrons can more easily move among the sites and 
their kinetic energy (which favors FM) is a stronger con- 
tribution to the total energy. Seen another way, the 
application of periodic boundary conditions to a cluster 
with open boundary conditions effectively raises the den- 
sity of the cluster's environment, since the cluster then 
appears to be surrounded by other clusters of the same 
density. 

A comparison between odd-A^ and even-A s clusters 
shows that clusters with an odd number of sites (which 
have integer spin for ±le~ away from half-filling) gener- 
ally have greater average spin relative to the minimum 
possible spin (zero for ±le _ ). This difference is not great, 
however, and their absolute average spin (e.g. in Fig. [22]) 
is comparable to that of the even- N s clusters, which have 
a minimum spin of 1/2 as opposed to 0. 

Cluster size is a third point of comparison, where we 
find that larger clusters usually have ground states with 
higher spin, and higher average spin overall. One should 
keep in mind, however, that larger clusters are able to 
have higher spin values just by virtue of having more 
sites (and total electrons). (Indeed, we find that smaller 
clusters have larger average spin relative to their maximal 
allowed spin.) The rise in average spin and the existence 
of higher spin ground states as cluster size increases is 
greater and more consistently true of electron-doped clus- 
ters. In this case the dependence of average spin on clus- 
ter size is particularly significant: we find a substantial 
percentage of maximally polarized clusters for all sizes 



(4-7) investigated, showing that the spin polarization in- 
duced by extra electrons persists to larger random sys- 
tems yielding large spins (up to S = 3). We also see that 
the polarization of larger clusters (6-7 sites) remains (and 
sometimes increases) when there are two electrons above 
half- filling. The average spin of hole-doped clusters shows 
a much weaker shift toward larger spin values with cluster 
size than the electron-doped case, which again highlights 
our central argument that electron-doping is very differ- 
ent from hole-doping. 

Fourth, we see that with increasing density there are 
usually fewer high-spin clusters in all categories except 
for clusters with one extra electron that have ty set by 
method 2 above (t > t). In this case the distribution with 
highest weight on large spins occurs at intermediate den- 
sity (p — in 2D, Tyrr in 3D), a result also seen in the 
ensembles of section \V\ below. This suggests that there 
exists an optimal density for finding high-spin states in 
doped semiconductors above half-filling. We generally 
expect low density to be most favorable for FM, since this 
corresponds to large U/t, and believe this is the reason 
why all but the aforementioned case show this behavior. 
In the exceptional case, when there is one extra electron 
and tij is set by model 2, the additional parameter t/t, 
will play a significant role, and the dependence of the 
pair (U/t, t/t) on the density could result in an optimal 
density for FM that is greater than zero. 

Lastly, the most striking trend we find is by compar- 
ing electron-doped and hole-doped clusters. When t = t 
the clusters with one extra electron have a spin distri- 
bution shifted to substantially higher spin values than 
those with one less electron (i.e. one hole). When is 
determined by our band calculation (i.e. > Uj), this ef- 
fect increases dramatically (particularly at intermediate 
density as mentioned earlier). This effect is expected, 
since in our model an extra electron hops with ampli- 
tudes tij while an extra hole hops with amplitudes ty. 
Recall that the motivation for the model comes from 
the special properties of the hydrogen atom which result 
in mobile electrons having spatially larger wavefunctions 
than mobile holes. These cluster results show that even 
in strongly disordered systems a Nagaoka-like ferromag- 
netism can emerge at least on the nanoscale, and one of 
the ideal conditions for this FM is an electron-doped sys- 
tem. Compared to those below half-filling, systems above 
half-filling also hold greater promise for spin polarization 
on longer length scales, since this would most likely arise 
from many aligned high-spin clusters. 



2. Large t case: t = U/2 

In model 3, the hopping iy is set to a constant C = 
U/2, a value near the maximum of t(r) (used in model 
2). This corresponds qualitatively to the case when the 
wavefunction on doubly-occupied sites is extended across 
the system (as if, for instance, the state had merged with 
conduction band states). One may access this regime ex- 
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FIG. 22: (Color online) Ground state average spin of 2D random clusters with fixed size and density, and open boundary 
conditions, as a function of electron-doping (negative = hole-doping). The lower half of plots are the result of setting fcy = ty, 
determined by the bandwidth of the lower Hubbard band. The upper half use t%j determined by the bandwidth of the upper 
Hubbard (D~) band. 



perimentally if the binding energy of the D~ state can be 
tuned (e.g. in many- valley semiconductors or by an ap- 
plied field). Here we focus on the case of 2D clusters with 
open boundary conditions, for which the raw spin distri- 
bution data comparing models 2 and 3 is shown in Fig. [T] 
Data for 3D clusters can be found in Appendix [Cl Two 
trends found in our discussion of models 1 and 2 above 
also appear in the i — U/2 results: odd-7V s clusters have 
greater spin polarization relative to their minimum spin, 
and spin polarization increases with cluster size. Unlike 
the results of method 2 with one electron (le), where 
the intermediate density was optimal for FM, the results 



of method 3 show spin polarization increasing with de- 
creasing density (as in method 1). This fits with our be- 
lief that the optimal density found when method 2 was 
used is due to the interplay of two density-dependent 
Hubbard model parameters (in method 3 there is only 
one, U/t, as in method 1). The electron-hole asym- 
metry found when t — U/2 is qualitatively similar to 
when tij — t(rij) (method 2), but with higher spin val- 
ues (for both electron- and hole-doped systems). This 
is expected in the single electron-doped case (le), since 
the second electrons are even more weakly bound, caus- 
ing correspondingly stronger spin polarization. In the 
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FIG. 23: (Color online) Ground state average spin relative to minimum spin of 2D random clusters with fixed size and density, 
and open boundary conditions, as a function of electron-doping (negative = hole-doping) . The lower half of plots are the result 
of setting Uj = Uj, determined by the bandwidth of the lower Hubbard band. The upper half use ty determined by the 
bandwidth of the upper Hubbard (D~) band. 



hole-doped case two aspects are particularly noteworthy. 
First, we find that the large t — U/2 results in clusters 
with higher spin than those of method 2, opposite to 
the trend seen in hole-doped bipartite lattices (cf. Figs.[S] 
and [T2")) . Second, the largest spin in the distribution sat- 
urates at a value of one below the maximal allowed spin, 
denoted S max (for instance, in 2D clusters with N s = 6, 
the spin distribution in nearly 100% S — 1.5). This be- 
havior is somewhat similar to the hole-doped triangular 
lattice (Fig. [T4|) , which has a partially-polarized ground 
state (with spin = S max — 1) which covers larger inter- 
vals of U/t as t/t is increased. In summary, large t results 



in almost 100% of (single) electron-doped clusters hav- 
ing ground state spin S max , and almost 100% of (single) 
hole-doped clusters having ground state spin S max — 1. 



V. CLUSTER ANALYSIS OF LARGE SAMPLES 

We now study the viability of ferromagnetism in a 
macroscopic sample. For this section, because of its rele- 
vance to hydrogenic n-doped semiconductors, we present 
results for and from the band calculation only. 
Our strategy will be to consider a large two- or three- 
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FIG. 24: (Color online) Percentage of magnetic clusters (spin 1 or greater) in an ensemble of 2D random clusters with fixed size 
and density, and open boundary conditions, as a function of electron-doping (negative = hole-doping). The lower half of plots 
are the result of setting ty — tij, determined by the bandwidth of the lower Hubbard band. The upper half use tij determined 
by the bandwidth of the upper Hubbard (D~) band. 



dimensional system of random sites and divide it into 
clusters that can be approximately treated as indepen- 
dent as far as the Hubbard part of the Hamiltonian is 
concerned. Choice of the number of carriers in each clus- 
ter involves long-range Coulomb forces and is treated in 
a classical approximation described later. We solve the 
clusters individually and then analyze the resulting dis- 
tribution of their ground state spins. The analysis of 
section IIVDI characterized random clusters with a fixed 
density; here the average density of a large system is 
fixed while the local density of individual clusters is free 
to vary. 



A. Decomposition into clusters 

We begin with a set of N sys randomly positioned points 
with some average density p where N sys is typically 
10,000 to 1,000,000. We then divide the points into ap- 
proximately isolated clusters, solve the cluster general- 
ized Hubbard Hamiltonian exactly, and consider their 
ground state statistics. We choose to divide the large 
set of points into clusters using a simple algorithm that 
proceeds as follows: 

1. Initially each point is a single cluster, and all points 
are "unused" . 
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TABLE I: Comparison of large t — U/2 and band calculation t distributions of ground state spin values for 2D random clusters 
with open boundary conditions. Table entries give the percentage of clusters with the ground state spin specified in the column 
header. Results are the ensemble average of many clusters with fixed size N s , density p, and doping = one electron (le) or hole 
(lh). Estimated error ±0.5%. 



2. Choose any unused point p, and find its nearest 
neighbor q. 

3. Merge the cluster containing p with the cluster con- 
taining q. 



4. Set point p to "used" status. 

5. Repeat at step 2 until no unused points 



remain. 



In this way we form the smallest clusters such that each 
point belongs to the same cluster as its nearest neighbor 
(i.e. the point most strongly coupled to it). Note also 
that the minimum cluster size is 2. The advantage of 
this "nearest-neighbor" method is that it always keeps 
nearest neighbor points in the same cluster, which is de- 
sirable from a perturbation theory standpoint. It does 
not, however, guarantee that the clusters include all the 
hopping amplitudes of the original system above some 
threshold. We show in Fig. I2"57 a) the decomposition of a 
2D system into clusters using the algorithm. A weakness 
of the nearest neighbor method is that it will form sepa- 
rate clusters of strongly-coupled pairs even when they 
are nearby other clusters, and it is clearly seen from 
Fig. [2S7 a) that some of the neglected bonds are stronger 
than other bonds that are kept. On the same set of ran- 
dom sites, the result of an alternate algorithm that keeps 
all hopping amplitudes greater than a certain threshold 
(chosen so that the size of the clusters is not too large) is 
shown in Fig. 125T b) . This technique removes the problem 
of isolating strongly coupled pairs/triangles from other 
nearby sites, but it has the disadvantage of being very 



sensitive to the threshold, adding another degree of ar- 
bitrariness. We find that both methods give reasonable 
decompositions into clusters, and the choice of algorithm 
not unique. In this work, we use the nearest-neighbor 
method outlined above, and leave a more detailed as- 
sessment and comparison of clustering methods for later 
work. 

We first determine, for fixed average densities p = , 
thtt and tItt, the distribution of cluster sizes which con- 

160 160 ' 

verges to the density-independent values shown in Table 
HIl By considering clusters with < 8 sites, which are 
within the reach of exact diagonalization techniques, we 
can account for over 97% of the sites. The remaining 
large clusters are converted into smaller clusters (< 8 
sites) by removing the smallest number of weakest links. 

We can estimate the local density, p\ oc , of an iVg-site 
d-dimensional cluster with sites at positions fj, i = 1...N 
from the formula: 



Plo 



N, 



in 3D 
in 2D 



(16) 



where R c i, the average radius of the cluster, is given by 



Rri — 
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i Ns 



(17) 



(18) 



For clusters of a given size N s and electron number N e , 
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FIG. 25: (Color online) Example of decomposing a 100- 
site system into clusters. Part (a) uses the nearest-neighbor 
method, and part (b) the threshold method (both described 
in the text). The blue lines link points in the same clusters 
(not all hopping links between the points are shown). 



the local density will also have some variation about its 
mean. We plot the local density distribution of clus- 
ters with 2-7 sites for normalized global average density 
p = 1.0 in Fig. [26j We find that clusters of larger size 
have a lower mean density, that is, at lower local densi- 
ties the process of following nearest neighbors links has 
greater probability of connecting together a larger num- 
ber of sites. The reason for this trend is due to the sup- 
pressed probability of finding a group of mutual nearest- 
neighbors at low densities. Let us consider the simplest 
case of two sites and compare the probability of finding 
a nearest neighbor at distance r corresponding to local 
density p\ oc — r~ d with the probability of finding a near- 
est neighbor at this distance that also has the original 
site as its nearest neighbor. We call such points "mutual 
nearest neighbors," and the differential probability distri- 
bution is found by multiplying the probability of finding 
a nearest neighbor by the probability that the second site 
does not have a NN closer than the first site. We thus 





Percentage of clusters 


N s 


2D 


3D 


2 


22.9 


20.9 


3 


28.2 


25.0 


4 


22.0 


20.6 


5 


13.7 


14.7 


6 


7.2 


8.6 


7 


3.5 


4.8 


8 


1.5 


2.7 


9 


0.6 


1.3 


10 


0.3 


0.6 


>10 


0.1 


0.8 



TABLE II: Distribution of cluster sizes in a large 2D or 3D 
system of random sites with a fixed average density. Clusters 
are formed from smallest sets of sites such that each site is in 
the same set as its nearest neighbor. Since this criterion does 
not depend on the value of the average density, this table is 
valid for all fixed average densities. 
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FIG. 26: (Color online) Individual density distributions for 2- 
to 8-site clusters in two and three dimensions when the aver- 
age density p = 1.0. Note that the majority of the weight falls 
above p, indicating that the clusters chosen are significantly 
more dense than the average. The inset shows the long tail 
of the 2-site cluster curve, indicating the existence of strong 
pairs. 
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FIG. 27: (Color online) Probability density for finding a NN 
compared to that of finding a mutual NN (a NN that also has 
the initial site as its NN) vs. local spatial density p\ oc . The 
short-dashed line (the difference) shows the probability that 
a nearest neighbor is not a mutual nearest neighbor, and thus 
will lead to a cluster of > 2 sites. The average density is set 
to unity. 
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FIG. 28: (Color online) Average ground state spin vs. local 
density of 2D 5-site clusters. The pertinent range of local 
densities is divided into bins, and bar heights indicate the 
average ground state spin of the 5-site clusters whose density 
falls within the corresponding density bin. This data is from 
p — yip clusters, but the behavior is typical (see text). 



define the differential probability of finding a mutual NN 
at distance r by PmutuaiNNM = Pnn(r) * (1 - P nn {r)), 
where 



Pnn(r) 
Pnn(r) 



1 — exp — 



n d/2 nr d 

r(| + i) 



2^/_ 2 

r(f) 



nr 



exp 
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(19) 



(20) 



The function p nn {r) is the probability of finding a site's 
nearest neighbor between r and r + dr, and P nn (r) — 
Jo" Pnn{r')dr' is the probability of finding a pair with 
length less than or equal to r. As shown in Fig. [27] for 
2D, at large p\ oc the distributions p nn (r) and p m utuaiNN 
approach one another, indicating that most nearest- 
neighbor links form mutual NN pairs. However, at lower 
pioc, due to the more rapid decrease of p m utuaiNN(?") at 
large r, the distributions separate and there is greater 
probability that a nearest neighbor will not be mutual, 
and thus lead to a larger (at least size 3) cluster. The 
peak in p mu tuaiNN near 3.0 coincides with the peak in 
the probability of 2D 2-site clusters in Fig. [23 as one ex- 
pects. Clusters of greater than 2 sites will have greater 
probability density at lower pi oc , closer to the peak in 
Pnn(r)- PmutuaiNN near 1.75 (also shown in Fig. [27|). The 
distribution PmutuaiNN^) decreases much more rapidly at 
large r (small pi oc ) than p nn (r) does, as shown in Fig. 1271 
We diagonalize all the cluster Hamiltonians individ- 
ually, and compile the resulting data to arrive at the 
distribution of spin values from the ensemble of clusters 
(obtained from many different large system realizations). 
In this case, there is substantial fluctuation in the local 
density of clusters. Even the mean local density for clus- 
ters of different size is different; only the average density 
of the entire system is fixed. The results show the same 



general trends as the clusters with fixed local density de- 
scribed earlier. For comparison, the average spin and per- 
centage of magnetic clusters in two an d th ree dimensions 
are presented in Appendix [C] and Ref . Il07l respectively. 

We also find that there is a weak correlation between 
local density (pi oc ) and average ground state spin (S). 
We observe quite generally that 2D clusters with one 
extra electron (N e = N s + 1) have a peak in (S) near 
Ploc ~ 0.015 while those with one hole (N e = N s — 1) 
have relatively smaller values of (S) that are less sensi- 
tive to changes in pi oc . Figure |2"51 shows this typical be- 
havior for 5-site clusters with p = and N e = N s ± 1. 
Similar qualitative behavior is found for other clusters 
sizes 4 < N s < 7 and from systems with p = , j|p , 
though (S) tends to be higher for larger size clusters. 
The location of the peak at pi oc 0.015 is important 
to our consideration of different large-system densities p, 
since the density-independent histogram of local density 
given in Fig. l2"fJl shows that clusters with p\ oc /p G [2,4] 
are most prevalent. In the case p = = 0.00625, 
Pioc = 0.15 corresponds to p\ oc /p — 2.4, whereas for 
P = TeTjo = 0.000625 and p = ^ = 0.01875 the cor- 
responding values of pi oc /p are 24 and 0.8 respectively. 
This suggests that the p — case will show the great- 
est overall magnetism, an inference that was seen in the 
fixed density clusters of section IIVDI and is supported 
by the further investigation below (see section IV Bp . 

So far we have focused on characterizing the ground 
state spin distribution for clusters with fixed size and 
electron number (but with varying densities). We now 
turn to the spin distribution of the large systems from 
which we have taken the clusters. Consider a large sys- 
tem with a fixed number of sites N sys and doping (fixed 
total electron number Nj; ot ). The system is partitioned 
into clusters of size N s — 2 — 7, which are approximated 
as being independent. It only remains to determine how 
the electrons will be distributed among the clusters - af- 
ter the number of electrons on each cluster is known, the 
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TABLE III: Average energy (in units of U » 0.945Ry*) re- 
quired to add (+1) or remove (-1) an electron from a half-filled 
cluster of size sites in a large 2D system with total average 
density p. We have used the t(r) and t(r) (with t > t) of our 
band calculation. 
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TABLE IV: Net average energy (in units of U) required to 
transfer an electron between a half-filled cluster of the size 
specified by a column to a half-filled of the size specified by 
the row. This data is for clusters in a large 2D system with 
total average density p — and with i(r) and t(r) set by 
our band calculation. 



clusters can be independently solved and their ground- 
state spin tabulated. We calculate the electron distribu- 
tion using three different methods, two of which ignore 
Coulomb interactions and one which takes them into ac- 
count using a classical approximation. In the following 
sections we consider only 2D systems, since our interest 
is primarily in 2D heterostructures and we can obtain 
better statistics in 2D than in 3D. 



B. Electron distribution without Coulomb 
interactions 

As a first attempt to find the distribution of electrons 
among the clusters, we ignore Coulomb interactions en- 
tirely and minimize the total energy, which is then just 
the sum of the cluster energies. To accomplish this, we 
must compute the ionization energy and electron affini- 
ties of (half-filled) clusters and minimize the system's to- 
tal energy to determine where the electrons will reside. 
We pursue this goal in two ways. 

a. Average energy method The first, more approx- 
imate, calculation finds the (ensemble) average energy 
required to remove or add an electron to a half-filled clus- 
ters of each considered size. These values are shown in 
Table IIIII The averages are over very broad distribu- 
tions, however, with standard deviations comparable to 
the mean value shown in Table IIIII This reveals a ma- 
jor shortcoming of this technique: it approximates broad 
energy distributions by their mean. Its advantages lie in 
the simplicity and speed of its calculation, and that it ap- 
plies to thermodynamically large systems. We continue 
the analysis, knowing that results are to be treated only 
as a first approximation. 

By subtracting pairs of the values in Table Hill we find 
the average net energy gained (or lost) when transferring 
an electron from one cluster to another, shown in Table 
IIVI The fact that all transfer energies are positive implies 
the stability of the electron configuration in which each 
cluster is exactly half-filled in this approximation. Note, 
however, that Coulomb interactions (which lower the en- 
ergy of two charged cluster system) may alter this picture 



substantially. Using the average affinities and ionization 
energies, we determine the optimal distribution of elec- 
trons among the clusters. Let x^ be the fraction of the 
total clusters that have n sites and charge q. In our cal- 
culation, n = 2... 7 and q G { — 1,0, +1} (clusters are 
allowed at most one electron or hole on them), so there 
are 18 variables in all. The optimal x^ are found by 
minimizing the total energy, E tot ({x%}), subject to con- 
straints. The energy is written: 

E t ° t ({x? l })= £ a n , q xl (21) 

where a n -\ is the energy required to add an electron to 
a n-site cluster, a n ^+\ is the energy required to remove an 
electron from a n-site cluster. Constraints on the problem 
are: 

• x n > for all n, q. 

• Yl q =-i x n ~ /"! where /„ is the fraction of total 
clusters with size n (found from Table [TTJ) . 

• Nl ot = J2 n q( n + where the sum ranges over 
all n = 2 ... 7 and q e { — 1, 0, +1} (note that n + q 
is the total number of electrons on n-site clusters 
with charge q). 

Since the energy and all constraints are linear, this is a 
linear programming problem and can be solved with stan- 
dard numerical routines. The minimization is carried out 
at a fixed doping, and results in the optimal placement 
the electrons in a thermodynamically large system. We 
determine the average spin per cluster of the entire sys- 
tem as a function of doping (summing the average spin 
of a cluster with n sites and charge q multiplied by x^) in 
Fig. [29] We also consider the percentage of clusters that 
have above-minimal spin (again using the average results 
for fixed-size clusters), classified above as magnetic clus- 
ters. From the plot of the percentage of magnetic clusters 
vs. doping in Fig. [30] we see that when a system is doped 
10-20% above half-filling, nearly half of the clusters have 
greater than minimal ground state spin. This suggests 
that at such doping some kind of percolation might be 
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FIG. 29: (Color online) Average spin per cluster as a function 
of filling (number of electrons per site; half-filled corresponds 
to 1.0), where the energy optimizing electron distribution is 
used at each filling. 
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FIG. 30: (Color online) Percentage of magnetic clusters (those 
having above minimal ground state spin) as a function of fill- 
ing (number of electrons per site; half-filled corresponds to 
1.0), where the energy optimizing electron distribution is used 
at each filling. 



possible that would induce magnetic order on a meso- 
scopic or even macroscopic length scale. 

b. Full cluster method A more straightforward way 
of calculating the optimal electron distribution in the ab- 
sence of Coulomb interactions is to consider an ensemble 
of large random systems. For each system, after sepa- 
rating the sites into clusters, we minimize the energy by 
repeatedly testing if the movement of an electron from 
one site to another lowers the total energy. Specifically, 
the algorithm we use is as follows: 

1. Initialize the system by placing electrons (if above 
half- filling) or holes (if below half-filling) on random 
clusters. 

2. Randomly choose two clusters i and j, and attempt 
to move an electron from i to j. If the resulting 
change in total system energy (just the sum of all 
cluster energies since there are no Coulomb inter- 
actions) is negative, accept the transfer. If not, do 
not make the transfer. 

3. Repeat the above step until the total energy con- 
verges. 



Knowing the electron distribution, the ground state spin 
of each cluster in the large system is then calculated. 
Finally, we compute the distribution of cluster ground 
state spin values, and average them over an ensemble of 
large systems. Figure [31] shows how the percentage of 
clusters with spin greater than or equal to a reference 
spin S re f. For S re f = 1, these percentages correspond 
to our earlier definition of "magnetic clusters" , and we 
compare in Fig. [32] the results of this section with those 
obtained earlier using the average energy method (section 
IV B al) . We see that the latter overestimates the number 
of high-spin clusters, particularly in the electron-doped 
case. The above analysis neglected the effect of long- 
range Coulomb interactions, which we consider next. 



C. Electron distribution including Coulomb 
interactions 

In low density insulating systems, Coulomb inter- 
actions between charged centers (or clusters) are not 
screened effectively and, because of their slow (1/r) fall- 
off, cannot be neglected*^ Therefore, a more accu- 
rate way to calculate the electron distribution is to in- 
clude Coulomb interactions. The approach described in 
this section accounts for the Coulomb interactions be- 
tween charged clusters in an approximate way. We be- 
gin in a fashion similar to the preceding analysis, solv- 
ing each cluster for range of total electrons near half- 
filling. We then determine the minimum energy elec- 
tron distribution by solving a generalized electron glass 
proble m 108 ! 109 ! 110 ! 111 which accounts for the differences 
in ground state energy of the clusters and the Coulomb 
energy between charged clusters, as described below. 

The generalized electron glass problem we solve con- 
sists of a set of two-dimensional clusters lying in a large 
two-dimensional space, indexed by i — 1 . . . N c i . Each 
cluster is treated as an effective site, and is assigned a 
position Ri (given by the average positions of all of its 
points), and a dimensionless charge qi. The charge nat- 
urally corresponds to the occupation of the cluster (rel- 
ative to half- filling) , and is restricted in our analysis to 
be +1, 0, or - 
Hamiltonian 



1. The problem is to minimize the classical 
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(22) 



where e is the dielectric constant, ry = \Rj, — Rj\, and cj)f' 
is the ground state energy of cluster i when it has charge 
qi. The first term gives the on-site (or, more accurately, 
on-cluster) energy contribution, and the second term sup- 
plies the Coulomb interaction between clusters. The min- 
imization is performed with respect to the variables qt 
which must obey the constraint Y2i 1i = A^* ot — N. 
where Nl ot is the total number of electrons in the N, 
site system. The details of the minimization are a gen- 
eralization of the procedure outlined by Baranovskii et 
al.r^ divided into three steps: 



sys, 
sys~ 
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FIG. 31: Percentage of clusters with total spin greater than 
the reference value S re f, specified in the key, as a function 
of filling (1.0 = half-filling). Plots correspond to densities 
P = TMo i lSo> an< l lfo as indicated in their titles. 
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FIG. 32: (Color online) Comparison of the percentage of mag- 
netic clusters (those with greater than minimal ground state 
spin) using the average energy method (AVG, thin lines) of 
section IV B al and the method using actual clusters without 
Coulomb interactions (NC, thick lines) of section lV B bl The 
style of the line (solid, dashed, or long-dashed) indicates the 
density. 



1. Initialize the {q{\ by starting them all equal to zero 
and randomly choosing clusters to add an electron 
to (if N* ot > N sys ) or remove an electron from (if 



Nl ot <N svs ) until Eli=Nl 



2. Calculate all single-cluster energies 



m 



and check that 



E 



Ef 



Ef 



< 



(23) 



(24) 



for all i, j such that qi > — 1, qj < 1, and i ^ j. The 
left hand side of the inequality is the cost of having 
the last-placed electron on site i, which should be 
less than the cost of placing an electron on site 
j. Otherwise, we can lower the system's energy 
(disregarding the Coulomb interaction for now) by 
moving an electron from i to j. In practice, we 
consider the pair i, j that for Eq. makes the left 
side maximal and right side minimal. If inequality 
(|24p is not satisfied we move an electron from i to 
j and repeat the step from the beginning. If the 
inequality is satisfied, we proceed to the next step. 
This is analogous to the /x-sub routine referred to 
by earlier work i 111 ! 112 

3. Calculate the energies Ei, and iterate through all 
pairs such that qi > —1, qj < 1, and i ^ j 

and check that each satisfies 



E 



- Ef - IK 



Ef 



> . (25) 



cr, 



If a pair is found that does not satisfy the 

inequality, we move an electron from cluster i to 
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cluster j and repeat the step (recalculate the Ei 
and check again). 

This process results in a set {qi} that is stable with 
respect to any single electron moving between clusters. 
Further conditions (and steps) could be added that would 
make the final distribution of electrons stable against any 
two electrons simultaneously moving between sites, but 
previous work on the electron glass problemiiS has shown 
that these additional constraints do not significantly af- 
fect the result. Therefore, we do not implement this ad- 
ditional step. 

Once we have determined the distribution of electrons 
among the many clusters of the large system, we com- 
pute the percentage of clusters with spin greater than a 
given reference spin S re f ■ This quantity is averaged over 
many random realizations of the large cluster system. 
The (ensemble-averaged) percentage of clusters with spin 
> S re f for S re f = h, 1, and | is shown in Fig. [33] for our 
standard densities p — — , , and y|p . For S re f = 1 , 
these percentages correspond to our earlier definition of 
"magnetic clusters" , and Fig. [M] compares the results 
with those of the previous section (jVB 0b[) which ne- 
glects Coulomb interactions but is otherwise identical to 
the calculation performed here. We see that Coulomb in- 
teractions slightly deplete the number of high-spin clus- 
ters, particularly in the electron-doped case. This also 
shows that even in the presence of long-range Coulomb 
interactions there is a sizable percentage of magnetic clus- 
ters at modest electron-doping. In order for the magnetic 
clusters to percolate in a strictly 2D system, they must 
account for 50% of the system, which is only attained at 
large filling factors (« 1.2 in the best case of p = j^q)- In 
3D, however, the percolation threshold is much lower, so 
a parallel calculation in a 3D or thick 2D system (which 
behaves as a 3D system on short length scales) may yield 
even more promising results. We also remark that as 
the impurity density is increased at fixed doping the av- 
erage number of magnetic clusters has a non-monotonic 
behavior. There is an optimal impurity density (nearest 
to p = in our data) that results in an on-average 
maximal number of magnetic clusters. Altogether, the 
presence of many high-spin clusters provides a necessary 
ingredient for ferromagnetism on macroscopic, or even 
mesoscopic, length scales. 
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FIG. 33: Percentage of clusters with total spin greater than 
the reference value S re f, specified in the key, as a function 
of filling (1.0 = half-filling) . Plots correspond to densities 
p = jm™, j4q , and as indicated in their titles. 



VI. CONCLUSIONS 



We have formulated a Hubbard model appropriate 
for doped semiconductors, which has an occupation- 
dependent hopping term and therefore intrinsic electron- 
hole asymmetry characteristic of the hydrogenic center. 
This generalized disordered Hubbard model is numeri- 
cally solved using exact diagonalization on 2D finite lat- 
tices, selected symmetric clusters, and completely ran- 
dom clusters in two and three dimensions. We summarize 
the results of each in turn. 
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FIG. 34: (Color online) Comparison of the percentage of mag- 
netic clusters (those with greater than minimal ground state 
spin) when Coulomb interactions are ignored (section lV B b[) 
or accounted for (section IV CI using a generalized Coulomb 
glass analysis. The plot shows, for densities p = , and 

percentages of the no-Coulomb (NC) case and electron 
glass (EG). 



Our results on finite (periodic) lattices, as well as 
selected clusters and distorted/randomized versions of 
them, lead us to several important conclusions. First, 
high-spin ground states generally occur at large U/t (low 
impurity density) . On a bipartite lattice one carrier away 
from half-filling, Nagaoka's theorem guarantees a maxi- 
mal spin state in the limit U/t — > oo. In the finite lattices 
that satisfy Nagaoka's theorem, we find maximal spin 
states at large but finite U/t. In clusters (with less sym- 
metry than a lattice), high-spin ground states are found 
to be quite sensitive to the cluster geometry, though they 
all exist at large U/t. Second, the properties of the hy- 
drogen atom give rise to a crucial difference between the 
electron-doping and hole-doping of hydrogenic systems. 
In lattices as well as clusters we see a greatly enhanced 
occurrence of spin-polarization in electron-doped (above 
half-filling) systems. In systems above half-filling we also 
find that increasing t/t can significantly increase the like- 
lihood of this nanoscale ferromagnetism. These results 
confirm our expectation that the greater the spatial ex- 
tent of a doubly-occupied site's wavefunction (relative to 
the wavefunction of a singly-occupied site), the more fa- 
vorable spin polarization becomes. Lastly, we remark on 
the resilience of the high-spin ground states. By perturb- 
ing a cluster geometry that has a high-spin ground state, 
we find that larger values of U/t and t/t make the state 
more stable to geometric and random fluctuations in the 
hopping amplitudes. An assessment of high-spin state 
robustness is relevant to situations in which sites are in- 
dividually positioned within some tolerance. Overall, we 
have identified a regime where nanoscale FM exists (with 
some robustness) in finite lattices and artificially-made 
clusters of hydrogenic centers. 

The analysis of ground state spin behavior in com- 
pletely random clusters reveals several of the same con- 
clusions we found for the selected symmetric clusters. 



Namely, we find that electron-doping and a larger t/t 
favor spin polarization in random clusters as well. The 
electron-hole asymmetry found in all of the random en- 
sembles implies that in real semiconductor systems there 
is a significant difference between doping above and be- 
low half-filling. Spin-polarization is much more preva- 
lent in systems above half-filling, an effect which we 
again emphasize as arising from the physical properties 
of the dopant atom. Unlike in the case of selected clus- 
ters, where ferromagnetism is generally more prevalent at 
larger U /t, we find that within the low-density range con- 
sidered (well below the metal-insulator transition), there 
is an optimal density for finding high-spin (random) clus- 
ters. This interesting observation is likely due to clusters 
breaking up into separate, effectively disconnected, pieces 
at very low densities, which hinders carrier movement 
and thereby the alignment of spins in the ground state. 

We also study the problem of distributing electrons 
onto the cluster components of a large system. Of par- 
ticular interest is the relatively small effect of Coulomb 
interactions (between charged clusters) on the electron 
distribution. Coulomb interactions slightly decrease the 
number of clusters with above-minimal spin. This ef- 
fect was unexpected since Coulomb interactions reduce 
the energy cost of charged clusters, which generally have 
higher ground state spin than un-charged clusters. A 
detailed look at differences in the electron distributions 
with and without Coulomb interactions may help to ex- 
plain the reason for the small observed effect, and is left 
for future work. 

Taking into account all our data on finite systems, we 
expect high spin clusters to be observable in systems with 
a low density (large Hubbard U/t) of centers and a small 
excess of electrons. The latter requirement is difficult to 
realize in 3D bulk systems, but could be met in doped 
quantum dots and 2D heterostructures. For example, 
doped quantum dots with dopant number Nd = 6—15 
and a small excess of electrons N e — Nd = 1 — 2 would be 
ideal systems for finding high-spin ground states. Also, 
in modulated structures with dopants in both quantum 
wells and barrier regions, regions of excess electrons can 
be achieved, unlike in true bulk doped semiconductors. 
We also note that the artificial cluster geometries stud- 
ied in section TlV Bl have real world applications through 
recently developed technology which allows precise place- 
ment of phosphorous donors in silicon^ 

The same regime (low density, electron-doping) is also 
the most likely region for the possible appearance of true 
macroscopic ferromagnetism, as our calculations on the 
cluster constituents of large systems (in section [Vj) re- 
veal. Obtaining a conclusive answer to this question nu- 
merically, however, requires going beyond the small sizes 
possible with exact diagonalization methods. A possible 
route is to use more approximate methods such as density 
matrix or perturbative renormalization group methods 
in combination with other numerical techniques. Even if 
true ferromagnetism on the macroscopic scale is absent, 
we have shown that there should be a significant asym- 
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metry between the magnetic response of systems with ex- 
cess electrons above the half-filled (uncompensated) case, 
and those with a deficit of electrons from the half-filled 
case, (i.e. traditional compensated): the former should 
have a larger susceptibility in the paramagnetic phase at 
low temperatures. This prediction can be experimentally 
tested by an experiment that uses gates to tune the elec- 
tron density in a 2D layer. Also, if ferromagnetism is 
attained on large enough length scales, it may show up 
as hysteresis in transport measurements due to magnetic 
domains. Clearly the temperature scales at which these 
ferromagnetic tendencies will manifest temselvcs will be 
much below the scales for diluted magnetic semiconduc- 
tors like Galium Manganese Arsenide. This is due to 
several factors - (i) the energy scale for shallow impuri- 
ties is low; (ii) the ferromagnetism occurs only for large 
U/t, i.e. low dopant and carrier densities where J ~ t 2 /U 
is very small; and (iii) Nagaoka ferromagnetism in a Hub- 
bard band is a much subtler effect involving two compet- 
ing terms, compared with ferromagnetism arising out of 
a double exchange mechanism. Nevertheless, the demon- 
stration of high spin states and possible ferromagnetism 
in semiconductors doped with so-called "non-magnetic" 
shallow impurities would suggest that magnetism in semi- 



conductors is not limited to semiconductors with transi- 
tion metal elements, but is possible in a wider range of 
semiconductor based materials. 



VII. AKNOWLEDGEMENTS 

This research was supported by NSF-MRSEC, Grant 
DMR-0213706 and DMR-0819860. 



APPENDIX A: CLUSTERS BUILT FROM 
SQUARES AND TRIANGLES (ONE HOPPING 
PARAMETER): DETAILED PHASE DIAGRAMS 

The following table shows the ground state (T=0) 
phase diagrams of "non ring-like" clusters 4, 6, 8, and 
9 of Fig. [T5] The fixed electron number is given in the 
upper-right corner of each plot, and only non-trivial dia- 
grams and their electron/hole pair are shown (i.e. if the 
1-electron diagram is non-trivial, the 1-hole diagram is 
shown for comparison). 
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APPENDIX B: CLUSTERS WITH TWO 
HOPPING PARAMETERS: DETAILED PHASE 
DIAGRAMS 

Here we present details of the ground state phase dia- 
grams for the clusters listed in Fig. |T7l For each cluster 
with 0-2 electrons away from half filling (in either direc- 
tion), we have computed phase diagrams for the param- 
eter range t 2 /h G [1,10], h/U G [0.01,0.5], and t/t = 1, 
5, and 10. Due to the large number of diagrams, we 
have divided them into three roughly defined categories: 
trivial, simple, and complex. Trivial diagrams have no 
phase transitions (changes in ground state spin) within 
the considered parameter range. Simple diagrams have 
an abrupt transition between a single pair of ground state 
spin values, with very little parameter space where there 
are intermediate spin values (effectively, these diagrams 



contain a single phase boundary). Complex diagrams 
have substantial parameter space where the ground state 
spin ranges three or more values, so there is a sizeable 
region of itermediate spin. 

Diagrams categorized as either simple or complex are 
listed in table [V] by cluster number (see Fig. [T7|) and 
charge Q = N s — N e (-Q is the doping relative to half 
filling). The absence of data (for any cluster number 
1-23 and Q = — 2, — 1, 0, +1, +2) implies that the phase 
diagram is trivial, and the constant ground state spin 
can be found from Fig. 1171 For each simple diagram, 
table [V] gives its two predominant spin values (min and 
max) and the equation of a line with approximates the 
phase boun dary. Full diagrams for these cases are given 
in Ref. 11071 . Diagrams categorized as complex are indi- 
cated by a "C" in table[V] and are shown in their entirety 
in Figs. EMSZ1 
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TABLE V: Description of all non-trivial phase diagrams for clusters in Fig. 1171 Cluster geometry is referenced by (see 
Fig. I17p . and Q is the clusters charge, t/t is the ratio of both pairs of kinetic parameters, i.e. t/t — ti/ti = te/ti. In each 
case, the minimal and maximal spin, Smin and Smax are given, as well as a rough approximation of the region of maximal spin. 
If there is a substantial region of the phase space under consideration (see text) where the spin lies between S m in and S m ax, 
then a "C" (for "complex" classification) is placed in the final column, and the diagram is shown in Figs. I35H37I If the phase 
diagram is given in the main text, regardless of its complexity, an "M" is placed in the final column. 
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FIG. 35: Ground state phase diagrams that contain substantial regions of intermediate ground state spin (spin between the 
minimum and maximum attained in the explored parameter space). 




FIG. 36: Ground state phase diagrams that contain substantial regions of intermediate ground state spin (spin between the 
minimum and maximum attained in the explored parameter space). 




FIG. 37: Ground state phase diagrams that contain substantial regions of intermediate ground state spin (spin between the 
minimum and maximum attained in the explored parameter space). 
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APPENDIX C: RANDOM CLUSTER SPIN DATA 
1. Spin distribution data 

Spin distributions, given by the percentage of clusters 
with each possible spin, are given for two-dimensional 
clusters with open and periodic boundary conditions in 
Table IVIII In this table we only consider models 1 and 2 
(see text, section llVD|) for determining tij - a comparison 
of models 2 and 3 is given separately. To facilitate com- 
parison between the two types of boundary conditions, 
the difference between the two cases (left and right sides 
of Table lVII|) is shown in Table IVll The analog ous r esults 
for 3-dimensional clusters are provided in Ref. Il07l Also 
included there is spin distribution data comparing mod- 
els 2 and 3 for 3D clusters (analagous to the 2D results 
of Pig.II). 



2. Average spin and percentage magnetic clusters 
of fixed density clusters 



This section presents complete results for the average 
spin and the percentage of magnetic clusters for fixed 
cluster size as a function of doping (zero doping = half- 
filled). Results for 2D and 3D clusters with open and 
periodic boundary conditions are shown in the figures 
below, as indicated in their titles. 

In particular, Fig. [35J shows the average spin of 2D 
clusters with open and periodic boundary conditions, re- 
spectively. To remove an even-odd effect, Fig. [39] shows 
the same average spin but relative to S m i n (i.e. 0.5 is 
subtracted from cases of odd electron number). Figure 
l40l shows the percentage of magnetic clusters (defined as 
those with greater than minimal ground state spin) as a 
function of doping for the different cluster sizes. Corre- 
sponding results for three-dimensional clusters with open 
and periodic boundary conditions are given in Ref. Il07l . 
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TABLE VI: Difference between distribution of ground state spin values for 2D random clusters with open and periodic boundary 
conditions. Values are obtained by subtracting the sets of data in Table IVTT1 below. Thus, entries show increase or decrease in 
percentage when switching from open to periodic boundary conditions. Estimated error ±0.7%. 
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TABLE VII: Distribution of ground state spin values for 2D random clusters with open boundary conditions (left) and periodic boundary conditions (right). Table 
entries give the percentage of clusters with the ground state spin specified in the column header. Results are the ensemble average of many clusters with fixed size N s , 
density p, and doping = one electron (le) or hole (lh). t > t indicates that t is set by our band calculation (to be compared with the case t = t). Estimated error 
±0.5%. 
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FIG. 38: Ground state average spin of 2D random clusters with fixed size and density, as a function of electron-doping (negative = hole-doping). Data for systems 
with open and periodic boundary conditions is shown in the left and right table respectively. The lower half of plots are the result of setting tij = ty, determined by 
the bandwidth of the lower Hubbard band. The upper half use tij determined by the bandwidth of the upper Hubbard (D~) band. 
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FIG. 40: Percentage of magnetic clusters (spin 1 or greater) in an ensemble of 2D random clusters with fixed size and density, as a function of electron-doping (negative 
= hole-doping). Data for systems with open and periodic boundary conditions is shown in the left and right table respectively. The lower half of plots are the result of 
setting tij = tij, determined by the bandwidth of the lower Hubbard band. The upper half use tij determined by the bandwidth of the upper Hubbard (D~) band. 
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3. Average spin and percentage magnetic clusters 
from fixed density large systems 

We next consider large systems with a fixed number 
of sites N sys (10,000 to 1,000,000) and doping (7V e tot to- 
tal electrons) . Each system is separately partitioned into 
clusters of size N s = 2 — 7, which are approximated as 
being independent, and then diagonalized. The resulting 
data, averaged over many (~ 50) large systems, gives the 
ensemble average distribution of clusters' ground state 



spin. The same general trends appear here as for the 
clusters with hxed local density. For comparison we show 
the average spin and percentage of magnetic clusters in 
two dimensions in Figs. 04] and 221 (for each cluster size 
separately) . We use the same plot format as for the clus- 
ters of fixed density, but show only the open boundary 
condition case. Note that although the cluster size is 
fixed, there is substantial fluctuation in the local density 
of clusters; only the average density of the entire syst em 
is fixed. Results for 3D clusters are given in Ref. Il07l 
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FIG. 41: Ground state average spin of 2D random clusters (open b.c.) obtained from large systems (N sys = 1 X 10 6 ) with fixed 
average density p, as a function of electron-doping (negative = hole-doping). The lower half of plots are the result of setting 
tij = tij, determined by the bandwidth of the lower Hubbard band. The upper half use tij determined by the bandwidth of 
the upper Hubbard (D~) band. 
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FIG. 42: Percentage of magnetic clusters (spin 1 or greater) in an ensemble of 2D random clusters (open b.c.) obtained from 
large systems (N sys = 1 x 10 6 ) with fixed average density p, plotted as a function of electron-doping (negative = hole-doping). 
The lower half of plots are the result of setting Uj = Uj, determined by the bandwidth of the lower Hubbard band. The upper 
half use tij determined by the bandwidth of the upper Hubbard (D~) band. 
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